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COMPUTER PROGRAMS FOR PLOTTING CURVES WITH 
VARIOUS DASHED- LINE SEQUENCES 

By Robert N. Desmarais and Robert M. Bennett 
Langley Research Center 

SUMMARY 

Two FORTRAN-callable subprograms have been written to draw a smooth curve 
through a set of input points as a solid line or as a general sequence of long and short 
dashes. Subroutine LINSEQ draws conventional curves whereas subroutine CONSEQ 
draws smooth closed curves (contours). The subprograms are based on an approximate 
calculation of the arc length along the curve and spline interpolation along the arc length. 
Options are provided for smoothing of the input data and for offsetting the plotted curve 
from the input data points. The method of calculation of the arc length and the generation 
of the line sequence are described. Usage descriptions of the main subprograms, sample 
calling programs illustrating the various features of the subprograms, and sample plots 
are given. The subroutines should be readily adaptable to almost any computer-driven 
incremental plotter. 


INTRODUCTION 

When plotting several curves on one graph, a unique sequence of long and short 
dashes is commonly used to identify each curve. When such curves are hand-drawn, a 
reasonable degree of skill and judgment is required by the draftsman to produce a result 
of good quality. For example, some effort should be made to see that the sequence is 
properly terminated at the end of the line and to maintain constant dash lengths along the 
arc of the curve. The modern digital computer has provided the capability of generating 
such a large volume of results that quality hand plotting and fairing of curves is prohibi- 
tively slow. Thus, automatic plotting hardware and software have been developed. Con- 
ventionally, however, only solid curves are generated by connecting closely spaced points 
with straight lines. Apparently, little work has been done in the past on automatic com- 
puter generation of dashed-line sequences analogous to hand plotting, partly because of 
the complication of maintaining a constant dash length along the arc of the curve. 

The purpose of this paper is to describe a set of computer subprograms (called 
LINSEQ and CONSEQ) to generate dashed-line sequences for an automatic plotter. Sub- 
program LINSEQ draws a conventional (or open) curve which can be a multiple -valued 



function of either coordinate X or Y. Subprogram CONSEQ draws a curve that closes 
smoothly with continuous deflections and slopes such as would be used for contour lines. 

An option is provided for the smoothing of the input data such as is sometimes desired 
when plotting experimental data points. An option is also provided for offsetting the 
curve from the input data so that heavy curves can be drawn with a fine pen to reduce 
smearing of ink or so that nearly coincident curves with different line sequences can be 
distinguished. 

The method of calculating the arc length, the method of generating the dashed-line 
sequence, and sample calling programs with input data and output plots are presented. 
Listings of each subprogram and usage descriptions of the principal subprograms are 
given. The spline -interpolation subprograms used for curve fairing can be used inde- 
pendently of LINSEQ and CONSEQ, and their usage descriptions are also given. The 
programs were written in FORTRAN for use on the Control Data series 6000 computer 
systems using the Langley Research Center versions of the RUN compiler and SCOPE 3.0 
operating system. However, they should be adaptable to other computer-driven incre- 
mental plotters without undue difficulty. 

This paper is intended to serve both as a user's guide and as a programer's manual. 
Sufficient information for the user is contained in the descriptions of subprograms LINSEQ 
and CONSEQ in appendix A and in the examples in appendix B. The detailed descriptions 
of the algorithms used, contained in the body of the text, and the program listings in appen- 
dix C can be skipped by those not interested in modifying the programs or converting them 
so they may be used with other equipment. 

METHOD 

In order to fair a dashed curve through a set of data points, it is necessary to have 
some means of interpolating the data to generate the plot vector file used by the incre- 
mental plotter. It is also necessary to incorporate into the interpolation algorithm some 
scheme for insuring that the dashes are equally spaced along the curve. Both of these 
objectives have been realized in the programs described herein by basing the interpola- 
tion on a pair of cubic-spline interpolation functions, with the approximate arc length from 
the first data point used as the independent variable and with X and Y as separate depen- 
dent variables. 

Two different curve-fairing subprograms are available. Subprogram LINSEQ fairs 
a set of N data points with a smooth curve passing through the points starting at ^Xj,Yjj 
and terminating at (Xn,Yn) and having zero curvature at the two ends. Subprogram 
CONSEQ constructs a closed curve starting and ending at (XjjYj^ with no discontinuities 
in slope or curvature. The only essential difference between the two subprograms is that 
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in LINSEQ the interpolation is performed by using natural cubic splines (i.e., splines for 
which the curvature vanishes at the two ends) and in CONSEQ the interpolation is per- 
formed by using periodic cubic splines (i.e., splines for which the slopes, deflections, and 
curvatures are matched at the two ends). Because LINSEQ and CONSEQ are so similar, 
only LINSEQ will be described in detail. 

Organization of Subprogram LINSEQ 

The curve-fairing subprogram LINSEQ is organized in such a manner that it pri- 
marily scales the input data, configures the basic line sequence, and calls several sub- 
programs to perform specific tasks. The tasks performed by these subprograms are 
listed below and will be described in detail in succeeding portions of this report: 

Subprogram Task 


PRETRP 

Implements preinterpolation option 

ARCALC 

Computes arc length 

SMOOTHC 

Smooths input data if smoothing option is selected 

SPISET 

Initializes spline interpolation 

ARCPLT 

Generates a single dash of the curve 

CHLSKYS 

Solves simultaneous equations 

SPIFUN 

Spline interpolation (ordinates) 

SPIDER 

Spline interpolation (derivatives) 

LIMPLT 

Commands pen motion if within user-prescribed limits 

CALPLT 

Causes pen motion with pen up or down (This subprogram is 
not part of the LINSEQ package; it is the installation library 
routine that interfaces the user program with the plotting 
devices.) 
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The logical flow of subprogram LINSEQ is illustrated by the flow chart below: 



The meanings of the tested parameters NN, SI, and NSP (further described in appendix A) 
are as follows: 

NN number of input points to be plotted (K NN = 0, it is assumed that the "no" 

branch was executed on a previous call.) 

SI , approximate interval between interpolated points on the plotted curve in 

inches (If SI is set to a negative value, a preinterpolation to 2+NN-l points 
is requested.) 

NSP number of iterated smoothing passes applied to the input data points 

Note that LIMSET and ARCSET are initialization entry points of subroutines LIMPLT 
and ARCPLT, respectively. Also SPISET initializes certain parameters subsequently used 
by the spline interpolation functions SPIFUN and SPIDER called by ARCPLT. 
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Arc- Length Calculation 

A set of data points with i = 1, 2, . . N, is given. A third array Sj, the 

array of arc lengths, has to be computed; Sj is set to zero and succeeding elements of 
the array are computed from the equation 

Si = Si_i + ASi (i = 2, 3, . . .,N) 

where ASi is the computed length of the arc of the curve connecting data point 

which is labeled (i-1) in the sketch below, to data point (Xy,Yi^, which is 
labeled (i) in the sketch below. This length ASj[ is approximated by the arc of a circle 
of radius R connecting the two points. The radius R is obtained by averaging the 
curvature of the left and right circumscribed circles with radii r' and r, respectively, 
as shown in the following sketch: 



Since any radius in this sketch may be infinite, it is more convenient to work with half- 
curvatures (h, h', and hav). Let h=— , h’=-^, and hav = — • Then, 

2r’ 2r' 2R 


hav - 


h’ + h 
2 


and 


ASi = 


sin 


-1 


(b’hav) 


‘av 
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Since it is considered desirable to underestimate the arc length ASi, the arcsine function 
is replaced by its two-term Taylor series; the resulting equation is 


ASi = b' 


1 + 


(^b’hav)*^ 

6 


In the sketch on the preceding page, let A/2 be the area of the triangle (i-1, i, i+1) and 

A A ' 

h = -4- and h’ = • 


let A'/2 be the area of the triangle (i-2, i-1, i). Then, 
where 


b'bc 


b”b'c 


A = (Xi - X,_i)(Yi+i - Yi_i) - (Xi^l - Xi.i)(Yi - Yi_i) 

K = - ^if 

and 

- J(Xl+l - Xi.lf + (Yi^i - Y,.j)2 

The terms A', b', and c' are computed the same as A, b, and c except that all 
subscripts are reduced by 1. Also, b" is computed the same as b except that sub- 
scripts are reduced by 2. 

This procedure for computing arc lengths to data points is implemented in sub- 
routine ARCALC called by LINSEQ and by CONSEQ. The parameter IFLAG is used by 
LINSEQ and CONSEQ to determine the manner in which the ends are treated. For 
LINSEQ, the curvature is set to zero at the end points and For 

CONSEQ, there are no ends because the data are assumed to define a continuous closed 
curve. 


Spline Interpolation 

After the array Si, with i = 1, 2, . . ., N, has been constructed, the set of input 
data points ^Xi,Yi), with i = 1, 2, . . ., N, is interpolated by a pair of parametric inter- 
polation functions 

X = X(S) 

and 

Y = Y(S) 
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where S, the independent variable, is the arc length along the curve. The interpolation 
functions used are cubic splines (e.g., refs. 1 and 2). A cubic spline is a set of piece- 
wise cubic polynomials joined at each data point, such that ordinate, slope, and second 
derivatives match where the polynomials are joined. Two additional conditions are 
required to define a cubic spline uniquely. In the spline subprogram called by LINSEQ, 
the additional conditions are obtained by setting the two end-point second derivatives to 
zero. Such a spline is called a natural cubic spline. In the spline subprogram called by 
CONSEQ, the additional conditions are obtained by causing the end-point first and second 
derivatives to match. Such a spline is called a periodic cubic spline (provided the ordi- 
nates also match as they do in this case). 

The spline interpolation functions are computed by using a compact scheme based 
on an overrelaxation solution of a system of linear equations. This algorithm (described 
in detail in ref. 1) furnishes the second derivative at each data point. Nonconvergence of 
the overrelaxation solution would result in discontinuities in the first derivatives. Since 
such discontinuities would be unacceptable, the two spline initialization subroutines 
SPISET and SPISETP use the computed second derivatives to compute left and right first 
derivatives at the data points. Then the averaged first derivatives are passed to the 
spline evaluation subprograms SPIFUN and SPIFUNP. Now nonconvergence causes dis- 
continuities in the second derivatives, but these discontinuities are so small that they can- 
not be detected in the faired curves generated by LINSEQ and CONSEQ. 

Six subprograms are used to implement the spline interpolations. Subprograms 
SPISET, SPIFUN, and SPIDER are called by LINSEQ. Subprograms SPISETP, SPIFUNP, 
and SPIDERP are called by CONSEQ. These subprograms can also be used, independently 
of LINSEQ and CONSEQ, as general-purpose interpolation routines. 

Dash Configuration and Offset 

Six user-furnished parameters are used to configure a sequence of long and short 
dashes. These parameters and their meanii^s are listed below: 

SI approximate interpolation interval in inches (The actual interpolation inter- 

val DS is computed by the program.) 

LI number of intervals of length DS in the gap between dashes 

L2 number of long dashes per cycle (A cycle consists of a single string of long 

dashes followed by a single string of short dashes.) 

L3 number of intervals of length DS in each long dash 
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L4 number of short dashes per cycle 

L5 number of intervals of length DS in each short dash 

The dashed-line sequences below illustrate the use of the five line-sequence parameters 
(LI, L2, L3, L4, L5). For each line sequence, the interpolation interval SI is approxi- 
mately 0.1 inch. 


(1, 2, 6, 3, 2) 

- - (3, 1, 1, 0, 0) 

(1, 1, 5, 2, 2) 

( 0 , 0 , 0 , 0 , 0 ) 

(1, 1, 4, 0, 0) 


Note that if only a single-length dash is to be used, L4 (not L2) should be set to 
zero. Setting LI to zero causes a solid line to be drawn. 

The actual interpolation interval along the curve DS is computed in such a manner 
that the curve from point to starts and ends with a series of long dashes 

if LINSEQ is used. If CONSEQ is used, DS is computed so that there are an integral num- 
ber of cycles starting at ^Xi,Yi^. 

The line-sequence parameter LI, L3, or L5, as appropriate, is passed to subroutine 
ARCPLT which then draws a single dash of the prescribed length. The line-sequence 
parameters L2 and L4 are used as loop delimiters in LINSEQ. 

The offset parameter H is used to draw parallel curves. For example, if the pen 
diameter in inches was D and one wanted to draw a curve of width 5D, one could call 
LINSEQ five times with H set to -2D, -D, 0, D, and 2D. Whenever multiple calls to 
LINSEQ are used with the same data, the parameter NN (number of data points) should be 
set to zero for all but the first call. This procedure causes the first part of LINSEQ to be 
bypassed. A positive value of H causes the curve to be offset to the right. Occasionally, 
it is desirable to set H to ±D/2 to distinguish curves that coincide over a significant frac- 
tion of their lengths. 


Use of Frame Limits and Plotter Conversion 

The pen-motion instructions issued by LINSEQ or CONSEQ are transmitted through 
a short subprogram called LIMPLT. Subprogram LIMPLT scans the pen-motion instruc- 
tion and transmits it only if the resulting pen motion would be within user-prescribed 
limits. These limits are prescribed by the parameters XLIM and YLIM in the calling 
sequence to LINSEQ. The parameters XLIM and YLIM are two word arrays (in the same 
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units as the input X- and Y-arrays) used to delimit the frame within which pen motion is 
permitted. They are transmitted, in scaled form (inches), from LINSEQ to LIMPLT by 
a call to the entry point LIMSET. 

The use of frame limits in this manner facilitates the plotting of a subinterval of 
the input data while retaining the influence of the points that are off-scale during the 
interpolation process. It can also guard somewhat against the pen hitting the mechanical 
stops on a mechanical plotter during the debug stages of a program. 

Subroutine LIMPLT is very short and is the only subroutine of the LINSEQ/CONSEQ 
package that contains a pen-motion instruction to the plotting device. This feature makes 
conversion of LINSEQ from one plotter to another very simple. The version of LIMPLT 
described herein contains the installation -dependent plotting instruction. 

CALL CALPLT(X,Y,IPLM) 

which means that the pen will be moved to the point (X,Y) (expressed in inches) with the 
pen lowered if IPLM = 2 or with the pen raised if IPLM = 3. (The usage of CALPLT 
is further described in appendix B.) Suppose that one wished to convert this instruction 
for plotting on a cathode-ray tube using the instructions of the graphics display library 
for the Xerox Data Systems Sigma series computers. For this device the pen-down 
instruction is 

CALL LINE(X,Y,MODE) 


and the pen-up instruction is 

CALL BEAM(X,Y,MODE) 


where 


MODE = 0 
MODE = 1 
MODE = 2 
MODE = 3 


means long absolute 
means long incremental 
means short absolute 
means short incremental 


Then the instruction 


CALL CALPLT(X,Y,IPLM) 
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should be replaced by the two instructions 

IF(IPLM.EQ.2) CALL LINE(X,Y,2) 

IF(IPLM.EQ.3) CALL BEAM(X,Y,2) 

Preinterpolation and Smoothing 

The preinterpolation and smoothing options provide the user with some control over 
the manner in which the curve fairing is performed. 

The purpose of preinterpolation is to compensate for the fact that LINSEQ uses the 
arc length as the interpolation abscissa although, in many cases, it is more desirable to 
use X as the interpolation abscissa. This option is selected by adjoining a minus sign to 
the interpolation-interval parameter SI. Preinterpolation is performed by adding to the 
X-array (x^, with i = 1, 2, . , ., n) additional N - 1 abscissas 
i = 2, 3, . . ., N, obtained by averaging adjacent abscissas — that is, 

Xi_i +Xi 

Xi-1/2 = (i = 2, 3, . . ., N) 

and then computing the corresponding ordinates spline interpolating with X as 

the independent variable. Example 3 of appendix B shows how the quality of the faired 
curve can be improved by preinterpolation. If the X-array is not monotonically increas- 
ing, preinterpolation will be attempted with Y as the independent variable. If neither the 
X-array nor the Y-array is monotonically increasing, preinterpolation will not be per- 
formed; thus, the program itself usually inhibits preinterpolation when it is not appropri- 
ate. However, preinterpolation is not appropriate and is not inhibited by the program if 
the X-array is monotonic but dY/dX is infinite or nearly infinite within the interpolation 
interval. This behavior is illustrated by Example 4 of appendix B. The preinterpolation 
option is not available in the closed-curve subprogram CONSEQ. 

The smoothing option is intended for use when the input data have random errors 
such as those resulting from a measurement process. The option is selected by assign- 
ing a value greater than zero to the parameter NSP (number of iterations) of subprogram 
LINSEQ. The value of NSP defines the number of iterations of the smoothing procedure 
in which a locally centered, five-point least-squares cubic is used as the smoothing for- 
mula for X and Y as functions of the arc length S. That is, at each point the 
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values are replaced by computed values 

Xi = *(Sl) 
and 

Yi = g(Si) 

where i = 1, 2, . . N. 

For example, the X-smoothing function = 3^0 ^ cubic 

polynomial with the coefficients ^aQ, aj, a 2 , and a 3 ^ chosen to minimize the sum 

i+2 2 

j=i-2 >- 

The weights w^, which serve to localize the cubic, are taken to be 

J 

^ 2 

fi * (Sj - %) 

The reference length I is equal to K AS where K is the iteration number of smooth- 

— Sjj - Si 

ing passes 1 to NSP and where AS = — represents the average value of the arc 

distance between consecutive data points. Note that if the length I is small, the cubic 
is strongly localized and there is very little smoothing. The Y-smoothing functions g(S) 
are evaluated in the same manner. At the two ends, the five points are taken to be the 
five closest points rather than the ith point and two others on each side. The smoothing 
option is available both with LINSEQ and CONSEQ. The amount of smoothing is controlled 
by the number of iterations NSP. If NSP = 1, there will be very little smoothing. The 
effect of letting NSP be very large is to cause the interpolated curve to resemble a para- 
metric cubic over its whole range. Example 5 of appendix B illustrates the effect of 
smoothing for some particularly rough input data. 

The preinterpolation and smoothing options partially cancel each other and are not 
normally both used for the same curve. 
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CONCLUDING REMARKS 


A set of computer subprograms (called LINSEQ and CONSEQ) to generate dashed- 
line sequences for an automatic plotter are described, and their usage is illustrated with 
examples. Subprogram LINSEQ draws a conventional (or open) curve which can be a 
multiple -valued function of either coordinate X or Y, whereas subprogram CONSEQ draws 
a curve that closes smoothly such as would be used for contour lines. The options pro- 
vided for smoothing of the input data and for offsetting the drawn curve from the input data 
are also described. In addition, the method of calculating the arc length, the method of 
generating the dashed-line sequence, and the sample calling programs are presented. The 
subprograms should be adaptable, without undue difficulty, to computer-driven plotters 
other than the one used. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., December 21, 1971. 
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APPENDIX A 


SUBPROGRAM DESCRIPTIONS 

This appendix contains usage descriptions of subprograms LINSEQ and CONSEQ, 
the plotting subprograms that are the subject of this paper. It also contains descriptions 
of the other subprograms called by LINSEQ and CONSEQ that are essentially independent 
subprograms suitable for other applications. (Since PRETRP, ARCALC, ARCPLT, and 
ARCPLTP are special-purpose subroutines that primarily communicate with LINSEQ and 
CONSEQ through labeled common, usage descriptions are not given.) In addition, since 
LINSEQ and CONSEQ are intended to be compatible with the NASA Langley Research 
Center graphic output library, a description of subprogram ASCALE has been included. 
ASCALE describes the use of the interleave factor, adjusted minimum, and scale factor 
used by LINSEQ and CONSEQ. The examples in appendix B also call CALCOMP, - 
LEROY/BALLPT and CALPLT from the graphic output library, so these descriptions are 
included as well. 
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APPENDIX A — Continued 
Subroutine LINSEQ 

Language: FORTRAN 

Purpose: To draw a solid or dashed smooth curve through a set of data points. Various 

dashed-line sequences with long and short dashes can be plotted. 

Use : CALL LINSEQ(XX,YY,NN,K,XLIM,YLIM,NSP,H,SI,L1,L2,L3,L4,L5) 

XX, YY The names of the arrays containing the input values of X and Y to be 

plotted. LINSEQ expects the adjusted minimums in locations 
K*NN+1 and the scale factors in locations K*(NN+1)+1 (normally, 
NN+1 and NN+2, respectively, for K = 1) of XX and YY as described 
in ASCALE. XX and YY must thus be dimensioned at least 
K*NN+K+1 in the calling program (normally, NN+2 for K = 1). 

NN The number of input points to be plotted. For successive calls with 

the same XX- and YY-arrays (such as for different values of H or 
of LI to L5), NN should be set to zero to bypass the arc length cal- 
culation and interpolation, thereby redrawing the previous curve 
with revised line parameters. 

K The interleave factor of XX and YY if mixed arrays. (Normally, 

K=l.) 

XLIM,YLIM Two word arrays containing the lower and upper limits on X and Y, 
respectively, in units of XX and YY. These parameters allow the 
user to designate a rectangular frame, outside of which pen motion 
is not permitted. 

NSP The number of iterated smoothing passes applied to the input data 

points. A weighted, locally centered, five-point least-squares cubic 
is used. Smoothing does not affect XX and YY. (Normally, 

NSP = 0.) 

H The distance the drawn curve is offset normal to the input data in 

inches, positive to the right as the curve is drawn. (Normally, 

H = 0.) 
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APPENDIX A — Continued 


SI The approximate interval between interpolated points on the plotted 

curve in inches (0.05 inch is a reasonable value). The actual inter- 
polation interval DS along the curve is computed from SI so that the 
curve starts and ends with a series of long dashes. If SI is set to a 
negative value, a preinterpolation of the input data to 2*NN-1 points 
is requested. Preinterpolation will be performed only if either the 
XX- or YY-array is monotonically increasing. 

LI The number of intervals of length DS in each gap between dashes. 

L2 The number of long dashes per cycle. 

L3 The number of intervals of length DS in each long dash. 

L4 The number of short dashes per cycle. 

L5 The number of intervals of length DS in each short dash. 

Note: If LI = 0, L2 through L5 will be ignored and a solid curve will be 

drawn. If L4 = 0, L5 will be ignored and all dashes will be the 
same length. With these exceptions, LI through L5 must be posi- 
tive integers. 

Restrictions : The subroutine contains five internal temporary storage arrays: S(lOl), 
X(lOl), Y(lOl), DX(lOl), and DY(lOl). These arrays are used to store certain plotting 
parameters. Thus, the value of NN is limited to 101 (51 if preinterpolation is used). 
These five arrays are in labeled common blocks ARCPAR, SCALEX, SCALEY, SLOPEX, 
and SLOPEY, respectively. If the user intends to process more than 101 input points, 
he should declare these blocks, with a larger dimension of at least NN, in a program or 
subprogram that will be loaded ahead of LINSEQ. It might also be noted that interpola- 
tion and computations are made for input XX and YY that are outside the frame speci- 
fied by XLIM and YLIM. If input points are extremely far outside the frame, computa- 
tional time will be unduly increased. 

Method : The approximate arcs joining consecutive data points are computed. Then the 
input is interpolated by using natural cubic splines with distances along the arc as 
abscissas and with X and Y as ordinates. Because neither X nor Y is a privileged axis, 
this subroutine can be used to draw nonsingle-valued curves such as spirals. A draw- 
back of this feature is that more data points are needed for acceptable fairing than would 
be required if conventional interpolation (i.e., with X as the abscissa) were used. 
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Storage : LINSEQ 436s 

ARCPLT 2038 

PRETRP 1618 

SPISET 3008 

SPIFUN 12?8 

SPIDER 1268 

SMOOTHC 2208 

CHLSKYS 2248 

ARCALC 304g 

LIMPLT 718 

LABELED COMMON^ 
(ARCPAR,SCALEX, 
SCALEY,SLOPEX, > 774g 

SLOPE Y, and 
C SPISET) 

TOTAL: , 3634g 


Subprograms used: CALPLT and SQRT in addition to above list 
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Subroutine SPISET/Functions SPIFUN and SPIDER 
Language : FORTRAN 

Purpose : SPIFUN is a function subprogram to perform interpolation by using a natural 
cubic spline as the interpolation function. SPIDER is a function subroutine to compute 
the first derivative of a natural cubic spline. SPISET initializes certain parameters 
used by SPIFUN and SPIDER. 

Use : CALL SPISET(N,X,Y,D) 

YP = SPIFUN(XP,N,X,Y,D) 

DYDX = SPIDER(XP,N,X,Y,D) 

N The number of words in the input X- and Y-arrays. 

X The array of input abscissas. 

Y The array of input ordinates. 

D The array of N points used to store the first derivatives dY/dX at the data 

points. These derivatives are computed by SPISET and are required by 
SPIFUN and SPIDER. 

XP The abscissa at which interpolation or derivative evaluation is to be performed. 

If XP is not within the interval X(l) to X(N), linear extrapolation is performed 
by using the deflection and slope at the appropriate end. 

Three additional seldom used parameters JMAX, IFRMS, and RMS are transmitted via 
labeled common in the following manner: COMMON/CSPISET/JMAX, IFRMS, RMS 
The parameters JMAX and IFRMS are defaulted by a DATA statement. 

JMAX The number of iterations used to compute the derivatives (Default = 20.) 

IFRMS If 1, the RMS is computed; othermse, the RMS is not computed 

(Default = 0.) 

RMS The root-mean-square value of the discontinuity in the second derivative 

at each data point. 
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Restrictions : The X-array must be monotonically increasing. The arrays X, Y, and D 
must be dimensioned N or larger in the calling program. SPISET must be called 
before the first call to SPIFUN or SPIDER. 

Method: A natural cubic spline is completely defined for the interval (Xi_i,Xi) between 

t ^ • \ / 

data points by the four quantities Yj, and Yj^. Subprogram SPIFUN or 

SPIDER tests XP to find which interval it is in and evaluates the ordinate or slope 
(respectively) of the cubic defined by the deflections and slopes at the two ends of the 
interval. Subprogram SPISET computes, and stores in D, the slopes required by 
SPIFUN and SPIDER. In order to save storage, the system of equations defining the 
array D is solved iteratively as in reference (a) of this subroutine. Since the D-array 
is only an approximation, there may be discontinuities in second derivatives at the data 
points. The root-mean-square value of this discontinuity can be tested by the user if 
desired. 

Accuracy : The iterative scheme used to compute second derivatives is absolutely stable 
and converges at a rate of slightly more than two bits per iteration. The nominal num- 
ber of iterations is 20. If the program is used for plotting, where high accuracy is not 
required, it is recommended that the user change the nominal number to 6. 

Reference : (a) Greville, T. N. E.: Spline Functions, Interpolation, and Numerical Quad- 
rature. Mathematical Methods for Digital Computers, Vol. n, Anthony 
Ralston and Herbert S. Wilf, eds., John Wiley & Sons, Inc., c.1967. 


pp. 

156-168. 

storage; SPISET 

3008 

SPIFUN 

1278 

SPIDER 

1268 

COMMON 

3s 

Subprograms used: 

None 


Other coding information : These subprograms use reentrant code. Thus, independent 
calls to SPISET/SPIFUN/SPIDER can be used in the same program provided the actual 
parameter lists are distinct. 
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Subroutine SMOOTHC /Subroutine SMOOTHP 

Language : FORTRAN 

Purpose : To replace a set of empirical ordinates by smoothed values. Typically these 
subroutines are employed prior to interpolating data when it is suspected that noise 
may be present. Subroutine SMOOTHP is used when Y is a periodic function of X. 

Use : CALL SMOOTHC (Npc,Y,T) 

CALL SMOOTHP(N,X,Y,T) 

N The number of points in the X- and Y-arrays. 

X A one-dimensional array of input abscissas. If SMOOTHP is called, X(N)-X(1) 
is the length of the period. 

Y A one-dimensional array of input ordinates. On return, Y contains smoothed 
ordinates. Note that SMOOTHP sets Y(N) to Y(l). 

T A 25-word array of temporary storage. On input, T(l) must contain a number 
designating the degree to which the smoothing formula is locally weighted 
(see Method). Usually, T(l) is set to 1. 

Restrictions : The array X must be monotonically increasing. 

Method : Each abscissa Y(I), where I is the point index, is replaced by a new value com- 
puted by passing a least-squares cubic through the five points (X(J),Y(J)), with 
J = 1-2, . . ., 1+2. For the ends, Y(l), Y(2), Y(N-l), and Y(N), the five closest points 
are used when SMOOTHC is called. 

The least-squares cubic is computed by using weights W(J) = l./(S**2+(X(I)-X(J))**2) 
where S = T(1)*(X(N)-X(1))/N. 

The user can select the amount of weighting by prescribing T(l) on input. In the limit 
as S — 0, the effect of W(I) is to inhibit smoothing completely. In the limit as S — 
the smoothing formula becomes the conventional five-point least-squares cubic 
described in reference (a) of these subroutines. 

Reference : (a) Hildebrand, F. B.: Introduction To Numerical Analysis. McGraw-Hill 
Book Co., Inc., 1956, pp. 295-302. 
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Storage: SMOOTHC 2208 

SMOOTHP 244s 

Subprogram used: CHLSKYS 
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Subroutine LIMPLT 


Language : FORTRAN 

Purpose : To move plotter pen to new location with pen up or pen down if new location is 
within prescribed frame limits, but not to move pen if new location is outside pre- 
scribed frame limits. 

Use: 


Calling sequence: CALL LIMPLT(X,Y,IPEN) 

X,Y The floating-point values of page coordinates for pen movement 

(inches). 

IPEN = 2 Move to point (X,Y) pen down if (X,Y) within frame limits. 

IPEN = 3 Move to point (X,Y) pen up if (X,Y) within frame limits. 

Entry point: CALL LIMSET(XLIM,YLIM,IPEN) 

XLIMjYLIM One -dimensional arrays of dimension two containing the X and Y 

minimum and maximum page -coordinate values (inches) in loca- 
tions one and two, respectively. 

IPEN = 0 For limits to be set with call to LIMSET. 

Restrictions : A call to LIMSET with IPEN set to zero must be made prior to the use of 
UMPLT. 

Method : Point (X,Y) is checked and if XLIM(l) = X = XLIM(2) and 

YLIM(l) i Y = YLIM(2) are both satisfied, the basic pen-movement subroutine 
(CALPLT) is called with the value of IPEN passed to LIMPLT. If these conditions are 
not satisfied, the pen will not be moved but will go pen up to the next point that satisfies 
the conditions. Recycling can be accomplished with a subsequent call to LIMSET. This 
method permits a user to plot only the portion of a curve within the frame. 

Accuracy : The limits of the frame are approached only as closely as the last point within 
the frame. 
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Subroutine CHLSKYS 


Language ; FORTRAN 

Purpose : To solve a set of systems of linear algebraic equations AX = B, where A is 
a square symmetric coefficient matrix and B is a rectangular matrix of right-hand- 
side column vectors. On return, X (the matrix of solution vectors) is stored in B. 

Use : CALL CHLSKYS(A,N,B,M,NX) 

A The square coefficient matrix. It is assumed symmetric by the program, so 

the user need only compute the diagonal and upper triangle of A. 

N The order of A and the number of rows of B. 

B The matrix of right-hand-side column vectors. On return, the solution vec- 

tors are stored in B. 

M The number of column vectors of B. 

NX The column size of A and B as given in the dimension statement of the 

calling program. 

Restrictions ; Reliable results can be guaranteed only if the matrix A is positive 
definite. The dimensions of A and B in the calling program should be A(NX,NN) 
where NN = N and B(NX,MM) where MM = M. 

Method : Cholesky's triangular decomposition is used (e.g., see ref. (a) of this subroutine). 
That is, A is expressed as 

A = 

where S is an upper triangular matrix with ones on the diagonal, D is the diagonal, 
and S is the transpose of S. Then, successively, 

B’ = S-1b 

B" = D"1b' 
and 

B”’ =S"^B” 
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are computed. B”' is the matrix of solution vectors X. The decomposition of A 
and the three steps above are performed concurrently, so the matrices S, D, and S 
are never actually available in storage. The input matrix A is destroyed by the 
subprogram. 

Accuracy : This algorithm can fail catastrophically if A is not positive definite. If it 
is suspected that jAijAjij can be greater than |A^Ajj| for some values of i and j, 
a different algorithm (employing pivoting) should be used. Loss of accuracy can also 
result from ill-conditioning. While there is no provision in the subprogram for 
condition-number computation, a useful estimate of the condition number is given by 
C1*C2 where Cl = A(l,l) before the call to CHLSKYS and C2 = A(N,N) after the call 
to CHLSKYS. Note that the loss of accuracy due to ill- conditioning is almost algorithm 
independent. 

Reference : (a) Wilkinson, J. H.: The Solution Of Ill-Conditioned Linear Equations. 

Mathematical Methods for Digital Computers, Vol. H, Anthony Ralston 
and Herbert S. Wilf, eds., John Wiley & Sons, Inc., c.1967, pp. 65-93. 

Storage : CHLSKYS 224g 

Subprograms used: None 
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Subroutine CONSEQ 


Language : FORTRAN 

Purpose : To draw a solid or dashed smooth closed curve (contour) through a set of data 
points. Various dashed-line sequences with long and short dashes can be plotted. 

Use : CALL CONSEQ(XX,YY,NN,KpCLIM,YUM,NSP,H,SI,Ll,L2,L3,L4,L5) 


XX, YY The names of the arrays containing the input values of X and Y to be 

plotted. CONSEQ expects the adjusted minimums in locations 
K*(NN+1)+1 and the scale factors in locations K*(NN+1)+1 (normally, 
NN+1 and NN+2, respectively, for K = 1) of XX and YY as described 
in ASCALE. XX and YY must thus be dimensioned at least 
K*(NN+1)+1 (normally NN+2 for K = 1) in the calUng program. 
CONSEQ closes the curve by taking the last point to be (XX(1),YY(1)). 
Thus (XX(1),YY(1)) must not be the same point as (XX(NN),YY(NN)). 

NN The number of input points to be plotted. For successive calls with 

the same XX- and YY -arrays (such as for different values of H or 
of LI to L5), NN should be set to zero to bypass the arc length cal- 
culation and interpolation, thereby redrawing the previous curve 
with revised line parameters. 

K The interleave factor of XX and YY if mixed arrays. (Normally, 

K = 1.) 

XLIMjYLIM Two word arrays containing the lower and upper limits on X and Y, 
respectively, in units of XX and YY. These parameters allow the 
user to designate a rectangular frame, outside of which pen motion 
is not permitted. 

The number of iterated smoothing passes applied to the input data 
points. A weighted, locally centered, five-point least-squares cubic 
is used. Smoothing does not affect XX and YY. (Normally, 

NSP = 0.) 

The distance the drawn curve is offset normal to the input data in 
inches, positive to the right as the curve is drawn. (Normally, 

H = 0.) 


NSP 


H 
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SI The approximate interval between interpolated points on the plotted 

curve in inches (0.05 inch is a reasonable value). The actual inter- 
polation interval DS is computed from SI so that the curve contains 
an integral number of cycles. 

LI The number of intervals of length DS in each gap between dashes. 

L2 The number of long dashes per cycle. 

L3 The number of intervals of length DS in each long dash. 

L4 The number of short dashes per cycle. 

L5 The number of intervals of length DS in each short dash. 

Note: If LI = 0, L2 through L5 will be ignored and a solid curve will be 

drawn. If L4 = 0, L5 will be ignored and all dashes will be the 
same length. With these exceptions, LI through L5 must be positive 
integers. 

Restrictions : The subroutine contains five internal temporary storage arrays: S(lOl), 
X(lOl), Y(lOl), DX(lOl), and DY(lOl). These arrays are used to store certain plotting 
parameters. The value of NN is limited to 100. These five arrays are in labeled 
common blocks ARCPAR, SCALEX, SCALEY, SLOPEX, and SLOPEY, respectively. If 
the user intends to process more than 100 input points, he should declare these blocks, 
with a larger dimension of at least NN+1, in a program or subprogram that will be 
loaded ahead of CONSEQ. It might also be noted that interpolation and computations 
are made for input XX and YY that are outside the frame specified by XLIM and YLIM. 

If input points are extremely far outside the frame, computational time will be unduly 
increased. 

Method : The approximate arcs joining consecutive data points are computed. Then the 
input is interpolated by using periodic cubic splines with distances along the arc as 
abscissas and with X and Y as ordinates. 

Storage : CONSEQ 
ARCPLTP 
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SPIFUNP 

SPIDERP 

SMOOTHP 

CHLSKYS 

ARCALC 

LIMPLT 

LABELED COMMON^ 
(ARCPAR,SCALEX, 
SCALEY,SLOPEX, > 
SLOPEY, and 
CSPISET) 

TOTAL 


1368 

1418 

2448 

2248 

3048 

718 


7748 


37048 


(Note that CHLSKYS, ARCALC, LIMPLT, and LABELED COMMON are also 
used by LINSEQ.) 

Subprograms used: CALPLT and SQRT in addition to above list 


27 



APPENDIX A — Continued 


Subroutine SPISETP/ Functions SPIFUNP and SPIDERP 

Language : FORTRAN 

Purpose : SPIFUNP is a function subprogram to perform interpolation by using a periodic 
cubic spline as the interpolation function. SPIDERP is a function subroutine to com- 
pute the first derivative of a periodic cubic spline. SPISETP initializes certain param- 
eters used by SPIFUNP and SPIDERP. 

Use : CALL SPISETP(Npc,Y,D) 

YP = SPIFUNP(XP,N,X,Y,D) 

DYDX = SPIDERP(XP,N,X,Y,D) 

N The number of words in the input X- and Y-arrays. 

X The array of input abscissas. X(N)-X(1) is the length of the period. 

Y The array of input ordinates. Y(N) need not be furnished, as it will be set to 

Y(l) by the program. 

D The array of N points used to store the first derivatives dY/dX at the data 

points. These derivatives are computed by SPISETP and are required by 
SPIFUNP and SPIDERP. 

XP The abscissa at which interpolation or derivative evaluation is to be performed. 
If XP is not within the interval X(l) to X(N), the argument of the cubic is 
made to lie in that interval by adding or subtracting the appropriate multiple 
of the period. 

Three additional seldom used parameters JMAX, IFRMS, and RMS are transmitted via 
labeled common in the following manner: COMMON/CSPISET/JMAX, IFRMS, RMS 
The parameters JMAX and IFRMS are defaulted by a DATA statement. 

JMAX The number of iterations used to compute the derivatives (Default = 20.) 

IFRMS If 1, the RMS is computed; otherwise, the RMS is not computed 

(Default = 0.) 
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RMS The root-mean-square value of the discontinuity in the second derivative 

at each data point. 

Restrictions : The X-array must be monotonically increasing. The arrays X, Y, and D 
must be dimensioned N or larger in the calling program. SPISETP must be called 
before the first call to SPIFUNP or SPIDERP. 

Method: A natural cubic spline is completely defined for the interval between 

data points by the four quantities Yj_i, Y{_j^, Yj, and Yj. Subprogram SPIFUNP or 
SPIDERP tests XP to find which interval it is in and evaluates the ordinate or slope 
(respectively) of the cubic defined by the deflections and slopes at the two ends of the 
interval. Subprogram SPISETP computes, and stores in D, the slopes required by 
SPIFUNP and SPIDERP. In order to save storage, the system of equations defining the 
array D is solved iteratively as in reference (a) of this subroutine. Since the D-array 
is only an approximation, there may be discontinuities in second derivatives at the data 
points. The root-mean-square value of this discontinuity can be tested by the user if 
desired. 

Accuracy : The iterative scheme used to compute second derivatives is absolutely stable 
and converges at a rate of slightly more than two bits per iteration. The nominal num- 
ber of iterations is 20. If the program is used for plotting, where high accuracy is not 
required, it is recommended that the user change the nominal number to 6. 

Reference : (a) Greville, T. N. E.: Spline Functions, Interpolation, and Numerical Quad- 
rature. Mathematical Methods for Digital Computers, Vol. n, Anthony 
Ralston and Herbert S. Wilf, eds., Jolin Wiley & Sons, Inc., c.1967, 
pp. 156-168. 


Storage: SPISETP 

520s 

SPIFUNP 

1368 

SPIDERP 

1418 

COMMON 

38 

Subprograms used: 

None 


Other coding information : These subprograms use reentrant code. Thus, independent 
calls to SPISETP/SPIFUNP/SPIDERP can be used in the same program provided the 
actual parameter lists are distinct. 
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Subroutine CALCOMP 


Language : FORTRAN 

Purpose : This is the normal mode processor. The necessary parameters and linkage 
are set up to output a tape for the CalComp Model 780/763 Electro -Mechanical Plotter. 


Use : CALL CALCOMP 

Restrictions : This call must be given before the first call to a plotting routine. 
Storage : CALCOMP 546?8 total for all subprograms used 

Subprograms used : PLOTSW, PLT763, BLCR, STRC, TAPWRI, ENCODl, STRCALL, 
CREATE F, BOUNDCK, TRUNCL, and LOCATE 


Other coding information : The following is a list of messages, the circumstances under 

which they will appear, and the action taken: 

NO PLOTTING DEVICE SPECIFIED This message is printed in the output file and 

the job is ended in subroutine PLOTSW. 

This condition occurs if there is no initial- 
ization CALL CALCOMP in the program 
prior to using subroutines which generate 
plotting output. 

PLOTTING COMMENCED This message is printed in the dayfile when 

the first pen movement is encountered as a 
result of a program call to the plotting 
subroutines. 

the last CALCOMP BLOCK These messages are printed in the dayfile 

ADDRESS WAS xxx when the plotting is completed, xxx is the 

DATA PLOTTED = yyyyyy value of the last block address on the 

CalComp plotter tape. This block address 
is at the end of the last valid data, yyyyyy 
is the approximate number of data points. 
Plotting is completed either as a result of a 
CALL CALPLT(x,y,999) or when the job is 
ended due to an error recognized by the 
CalComp subroutines. 
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Subroutine CALPLT 


Language : FORTRAN 

Purpose : To move the plotter pen to a new location with pen up or down and to signal the 
end of a job segment by incrementing the block address number. 

Use : CALL CALPLT(X,Y,IPEN) 

X,Y The floating-point values for pen movement. 

IPEN = 2 Pen down 

IPEN = 3 Pen up 

Negative IPEN will assign (X = 0, Y = 0) as the location of the pen after moving the 
X,Y (creating a new reference point) and will increase the block number by one. 

(The block number is the number that appears in the display at the top of the tape 
drive on the plotter and identifies the portion of the output tape that is being plotted. 
The block address 001 is written automatically as a result of the initialization pro- 
cessor call.) Each block address generally implies a separate page or plot. 

IPEN = 999 Writes a terminating block address of 999 for peripheral handling of the 
plotter tape and all further processing is skipped. X and Y may be 
any values since they are ignored. 

Restrictions : All X- and Y- coordinates must be expressed as floating-point values in 
inches (actual page dimensions) in deflection from the origin. 

(A CALL TO CALPLT WITH EITHER NEGATIVE IPEN (USUALLY -3) OR A TERMI- 
NATING BLOCK ADDRESS (IPEN = 999) MUST BE GIVEN AS THE LAST PLOTTING 
INSTRUCTION BEFORE ENDING A PROGRAM WHICH USES ANY OF THE PLOTTER 
SUBROUTINES: THIS IS TO BE SURE THAT ALL PLOTTER INSTRUCTIONS ARE 
WRITTEN ON THE PLOTTER TAPE.) 

Method : The main subroutine in the CalComp software package is the CALPLT sub- 
routine. All other special-purpose subroutines eventually call CALPLT either directly 
or indirectly. Subroutine CALPLT moves the pen in a straight line between the present 
pen position and another pen location to which the programer wishes the pen to be 
moved. 
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In order to cause such instructions to be written, the programer specifies the coordi- 
nates of the point to which the pen is to be moved and whether the pen is to be moved 
in a raised or lowered position. This movement is accomplished by the FORTRAN 
instruction 

CALL CALPLT(X,Y,IPEN) 

Also, the subroutine provides "sequence numbers” on the tape, making it possible to 
afford identification of job segments. The block address 001 is written on the first 
call to CALPLT. Thereafter, if the programer defines a new origin or wishes to 
divide the job into several segments, he need only set the argument IPEN negative. 
The CALPLT routine then moves the pen to (X,Y); stores this location as (0,0), that 
is, a new origin; and increases the block address by one. 

Storage : CALPLT 2538 

Subprograms used: PLOTSW, STRCALL, and LCXJATE 
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Subroutine ASCALE 


Language : FORTRAN 

Purpose : To compute a scaling factor for an array of numbers to be plotted over a cer- 
tain area and find the minimum data value within the array. 

Use: CALL ASCALE(ARRAY,S,N,K,DV) 


ARRAY The name of the array containing the floating-point values to be scaled. 

S The length (floating-point value in inches) over which the data are to be 

plotted (usually the length of one of the axes). 

N The number of data values in ARRAY from which points are to be plotted 

in accordance with K. 

K The interleave factor which specifies the sequence in which data are 

stored. 

K = 1 indicates that values are stored sequentially. 

K = 2 indicates that values are stored in every other location in the array 
and so forth. 

DV The number of divisions per inch of the plotting paper to be used. (DV 

should be 10.0, 20.0, 25.0, or 25.4.) 

Restrictions : The array must be dimensioned to include storage space for two extra 
elements per interleave factor. For example: N = 100, K = 1, DIMENSION ARRAY 
(102); N=75, K = 3, DIMENSION ARRAY (231). 

Method : This routine scans the elements in the array to find the minimum and maximum. 
ASCALE computes an adjusted minimum (origin value) and stores it in ARRAY((N*K)+1) 
and computes a scale factor and stores it in ARRAY((N*K)+1+K). The scale factor will 
be either 2, 4, 5, or 10 with a decimal exponent. The data in the array may be scaled 
to floating-point values in inches by using a formula similar to the following: 

SV = (AE-MV)/SF, where SV is the scaled value, AE is the present value of the 
array element, MV is either the minimum value or the value desired at the origin, 
and SF is the scale factor computed by the subroutine. 
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Storage : ASCALE 262g 
Subprograms used : ALOG, ALOGIO 

Other coding information : Example: DIMENSION ORD(102),ABS(204) 

CALL ASCALE(ORD,10.,100,i,10.) 
CALL ASCALE(ABS,15.,100,2,10.) 
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Subroutine LEROY/BALLPT 


Language : FORTRAN 

Purpose : The parameters necessary to accommodate plotting with the liquid ink pen are 
set up by CALL LEROY. Once set, this mode will remain in effect as long as 
CALCOMP is in use or until a call to BALLPT is given. 

The parameters for plotting with the ballpoint pen are reset by CALL BALLPT. This 
mode is automatically in effect with CALCOMP unless there has been a call to LEROY. 

Use : CALL BALLPT 
CALL LEROY 

Restrictions : The CALL LEROY should be used only with CALCOMP. In addition to 
reducing the speed of the plotter for all plotting movements, the number of plot vectors 
in any annotation is considerably increased. 

The CALL LEROY must be made prior to any plotting calls, but after the CALL 
CALCOMP. 

Storage : LEROY/BALLPT 250 


Subprograms used: None 



APPENDIX B 


SAMPLE CALLING PROGRAMS AND THEIR RESULTS 


This appendix contains several examples illustrating the use of subroutine LINSEQ. 
For each example, the adjusted minimum was set to 0, the scale factor was set to 0.125, 
and the interpolation interval SI was set to 0.16 inch, except when otherwise noted. The 
curves were reduced to approximately one-half of their original size for insertion into 
this appendix, so the interpolation interval on the figures in this appendix is 0.08 inch. 

For each example, the five line-sequence parameters are listed in the form 
L = (L1,L2,L3,L4,L5). 

In order to keep the calling sequences simple, only LINSEQ and CONSEQ have been 
used to construct the plots. For this reason, most of the plots do not include the usual 
additional information, such as axes, scales, and legends. The sample calling programs 
also contain four installation-dependent calls that do not actually cause pen motion. They 
are as follows: 


CALL CALCOMP 
CALL LEROY 
CALL CALPLT(0.,0.,-3) 
CALL CALPLT(0.,0.,999) 


Device initialization 

Slows down device to improve quality 

Sets origin on paper 

Writes terminating block address on plotting tape 


Examples 3,4, and 5 also contain a call to a Langley subroutine LINPLT used to put 
plotting symbols at the input data points. The parameter list is (X,Y,N,K,- 1,13,3,0) 
where X, Y, N, and K have the same meaning as in LINSEQ; -1 means don't connect the 
points; 13 means use the 13th plotting symbol (diamond with cross inside); 3 means use a 
large symbol; and the last parameter 0 is unused here. 

Note that for the FORTRAN compiler used, the symbol $ separates distinct 
FORTRAN statements on the same coding line. 


Example 1 — Plotting Empirical Data 

Example 1 illustrates the use of LINSEQ to plot empirical data and shows that the 
input data need not represent a single-valued function. A letter D was drawn on graph 
paper and digitized at 15 points as shown in the figure and table on the following page. 
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X 

Y 

0.60 

0.65 

.15 

.30 

-.25 

-.35 

-.50 

-.65 

-.75 

-.60 

-.75 

-.50 

-.45 

o 

1 

-.05 

-.70 

.45 

-.60 

.60 

-.15 

.35 

.35 

-.10 

.70 

-.50 

.65 

-.60 

.15 

-.40 

.05 



Result .- The information in the above table was read from cards and plotted by 
using the line-sequence parameters L = (1,1, 3, 0,0). The resulting plot is 
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Listing of program and data .- The listing of the program and data for Example 1 is 
as follows: 


DIMENSION X( 101) ,Y( ton ,XLIM(2) ,YLIM(2J 
* example 1. PLOTTING EMPIRICAL DATA 
A0JMIN=0. i SCAFAC=.12? 

CALL CALCOMP $ CALL LEROY $ CALL C AL PL T ( C . 1 0 . 3 ) 

SI=.l6 t, K=1 
READ 2, N 

X(N+1)=Y{N+1)=A0JMIN i X(N+2)=Y(N+2)=SCAFAC 
XLI M( 1)=YLI M( 1)=-1.2 $ XLI M( 2)=YL IM( 2 ) = 1.2 

DO I 1=1, N 

1 READ 3, X( I ) , Y( I ) 

CALL LINSEO(X,Y,N,K,XLIM,YLIM,0, 0,51,1,1,3,0,0) 

CALL CALPLTIO. ,0.,9PP) 

2 FORMAT) I*-) 

3 FORMAT! 2F5. 2 I 
END 

1 '. 

.60 .65 

.15 .30 

-.25 -.35 
-.50 -.65 
-.75 -.60 
-.?5 -.50 
-.45 -.40 
-.05 -.70 
.45 -.60 
.60 -.15 
.35 .35 

-.10 .70 

-.50 .65 

-.60 .15 

-.40 .05 


Example 2 — Various Features of LINSEQ 

Example 2 illustrates the following aspects of the use of LINSEQ: 

(1) The use of LINSEQ to distinguish several different curves on the same plot. 

(2) The use of a user-written subroutine (GENSEQ here) to generate a "standard" 
set of line-sequence parameters. (This subroutine could be more complicated or simpler 
than GENSEQ.) 

(3) The use of two-point calls to draw axes and the identification key for the figure. 
(Note that the lettering Ui(X) through U 7 (X) was typed after the plots were 
generated.) 

(4) The use of user-declared labeled common blocks to increase the capacity of 
LINSEQ to more than the nominal value of 101 points. 
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APPENDIX B — Continued 


(5) The use of frame limits to limit the plotted curves to a predeclared frame. 
Result.- The data plotted are Chebyshev polynomials of the second kind 




sinl^L cos~^Xj 
2 sin (cos “lx) 


(L = 2, 3, . . ., 8) 


digitized at 23L - 43 points (i.e., from 3 to 141 points depending on degree). The frame 
limit parameter YLIM was used to limit the curves to Y = ±1. The resulting plot is 
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APPENDIX B - Continued 


Listing of program.- The program listing for Example 2 is as follows: 


DIMENSION X{ 143) ,Y( 143) tXLTM(2) ,YLIM{2) 

* EXAMPLE 2. VARIOUS FEATURES OF LINSEQ 
DIMENSION U(4),V(4) 

COMMGN/ARCPAR/ARCPAR(141)/SCALEX/SCALEX( 141 )/SCALEY/SCAL£Y( 141 ) 
COMMON/ SLOPFX/SLOPEXt 14 1 ) /SLOPE Y/SL OPEY ( 141 ) 

ADJMIN=0. $ SCAFAC=.125 

CALL CALCOMP $ CALL LEROY $ CALL CALPLT (0. ,0. ,-3) 

SI=.08 i K=1 

XLIM(l)=YLIM(l)=-1.0n001 $ XL I M( 2 ) =YL IMl 2 ) = 1 . 00001 

U(3)=V(3)=A0JMIN $ U(4)=V(4)=SCAFAC 

U(l)=-1. S U(2)=l. i V(1)=V<2)=0 
CALL LINSEO(U,V,2,K,XLIM,YLIM,0,0,SIfO,0,0,0,0) 

CALL LINSEQ(V,Ut2,K,XLIM,YLIM,0,0,SI»0,0,0,0,0) 

PI=3. 14159265358970 $ EPS=1.E-12 

DO 2 L=2,8 * N=23’«'L-43 $ DO 1 1 = 1, N 

T = PI4:{N-I)/(N-1.)+EPS $ X(I)=COS(T) 

1 Y(I )=.5*SIN(L*T)/SIN(T) $ CALL GENSEO ( L , L 1 , L 2 , L3 , L4 , L5 ) 

X (N + 1 )=Y (N + 1 ) = A0JMIN $ X(N4-2)=Y(N+2) = SCAFAC 

2 call LINSEQ(X,Y,N,K,XLIM,YLIM,0,0,SI,L1,L2,L3,L4,L5) 

VLIM(1)=-2. $ U(l)=-.^ $ U(2)=.4 $ V(1)=-1. $ DO 6 L=2,3 

CALL GENSEQ(L,L1,L2,L3,L4,L5) t V ( 1 ) = V ( 2 ) =V ( 1 ) - .05 

f CALL LINSE0(U,V,2,K,XLIM,YLIM,0,0,SI,L1,L2,L3,L4,L5) 

CALL CALPLT(C.,0.,999) 

END 


SUBROUTINE GE NSEO ( L , L 1 , L 2 , L3 ,L4 , L5 ) 


GENERATES A 
Ll=2 t. 

1 L1=0 $ 

2 L2=l S 

3 L3=16 $ 

RETURN 
END 


SO-CALLEO STANDARD LINE SEQUENCE FOR USE BY LINSEQ 

L2=L3=L4=L5=0 $ IF(L-2)1,2,3 

RETURN 

L3=6 $ RETURN 

L5=4 i IF(L.GT.5)L3=12 $ L2=L/3 * L4=L+1-3*L2 

SUBROUTINE GENSEO 


Example 3 — Preinterpolation 
Example 3 illustrates the use of the preinterpolation option. 

Result . - The function Y = sin(2irX) was generated at 21 equispaced abscissas on 
the interval (-1,1). These abscissas were plotted by using symbols. Then LINSEQ was 
called with line-sequence parameters L = (1,1, 2 ,0,0). The resulting curve that looks 
like a string of bowling pins (see short dashes) is much worse than one would expect from 
10 data points per period. LINSEQ was called again with line-sequence parameters 
L = (1,1, 6, 1,2) and with the preinterpolation option set (minus sign on interpolation inter- 
val SI). The curve that results from using preinterpolatLon (see long and short dashes) is 
very close to a true sine wave. These two curves are presented on the opposite page. 
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Listing of program .- The program listing for Example 3 is as follows 


DIMENSION X(23),Y(23» ,XLIM(2)fYLIM(2) 

FXAMPLE 3, PREINTERPQLATION 
A0JMIN=0. $ SCAFAC=.125 

CALL CALCOMP $ CALL LEROY * CALL CALPLT ( 0 . , 0 . ,-3 ) 
SI=.l6 $ K=1 $ N=2l 

X(N+1)=Y(N+1)=ADJMIN * X ( N+2 ) =Y ( N+2 » = SCAF AC 
XLIMI 1)=YLIM( 1)=-1.2 $ XLIM(2»=YLIM(2)=1.2 

PI=3. 14159265358979 
00 1 1=1, N $ X( n=-l.>, !♦( I-U 
1 Y(I ) = SIN(2.*PI*Xn ) J 

CALL LINPLT(X,Y,N,K, -1,13, 3,0) 

CALL LINSE0(X,Y,N,K,XLIM,YLIM,0,0,SI , 1, 1,2,0,0) 

CALL LINSEQ(X,Y,N,K,XLIM,YLIM,0,0,-SI,1,I,6,1,2) 

CALL CALPLT(0.,0.,999) 

END 



APPENDIX B — Continued 


Example 4 — Incorrect Use of Preinterpolation 

Example 4 shows how the preinterpolation option can be used incorrectly. 

Result .- A semicircle was generated at nine points and plotted as symbols. LINSEQ 
was called with line-sequence parameters L = (1,1, 2 ,0,0). The result is the short-dashed 
curve. Then LINSEQ was called again with L = (1,1,6, 1,2) and the preinterpolation 
option set. The resulting curve (long and short dashes) is much worse than the result 
without preinterpolation because dY/dY is infinite at the two ends of the interval. Pre- 
interpolation should not be used if Y or dY/dY is infinite or extremely large over the 
interpolation interval. These two curves are presented below: 



Listing of program.- The program listing for Example 4 is as follows: 


DIMENSION XI231,Y(23) fXLlM(2l,YLIM(2) 

* EXAMPLE 4. INCORRECT USE OF PREINTERPOLATION 
ADJMIN=0. $ SCAFAC=.125 

CALL CALCOMP $ CALL LEROY $ CALL CALPLT (0. ,0. i-3) 
SI=.16 $ K=1 $ N=9 

X(N+1)=Y(N4-1)=ADJMIN % X(N+2)=Y(N+2) = SCAFAC 

XLI M( 1)=YLI M( l)=-l .2 $ XLIM(2)=YLIM(2)=1.2 

91=3.14159265358979 

DO 1 I=lfN $ T=PI*( N-I l/IN-l. ) $ - X( I)=C0S( Tl 

1 Y( I ) = SIN(T) 

CALL LINPLT(X,YtN,K, -1,13, 3,0) 

CALL LINSE0(X,Y,N,K,XLIM,YLIM,0,0,SI,1 ,1,2,0,0) 

CALL LINSE0(X,Y,N,K,XLIM,YLIM,0,0,-Sl , 1, 1,6,1 ,2) 

CALL CALPLT(0.,0.,999) 

END 
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Example 5 — Smoothing 

Example 5 illustrates the use of the smoothing option. 

Result .- Twenty -two points from an empirical curve were read and plotted as 
symbols. LINSEQ was called three times with parameters set as shown below: 


Solid curve No smoothing (NSP = 0); L = (0,0, 0,0,0) 

Short dashes One smoothing pass (NSP =1); L = (1,1,3,0,0) 

Long and short dashes Eight smoothing passes (NSP = 8) ; L = (1,1, 6, 1,2) 


Note that the curve with 8 smoothing passes (NSP = 8) closely resembles a cubic polyno- 
mial in the resulting plot presented, below: 



43 



APPENDIX B — Continued 


Listing of program and data.- The program listing for Example 5 is as follows: 


DIMENSION X( 10l),Y( lOl) ♦XLIM(2» ,YLIM(2) 

* EXAMPLE 5, SMOOTHING 

CALL CALCOMP $ CALL LEPOY i CALL CALPLT (0. ,0. »-3) 
SI=.16 $ K=1 

READ 2, N 

XLIM(L>=YLIM(l)=-20. $ XLI M( 2 ) =YL IM ( 2 »=20 . 

READ 3, X(U,Y(n $ XMIN=XMAX=X( 1) t YMIN=YMAX=Vm 
DO 1 1 = 2, N 5 READ 3, X(I),Ym 

IF( X( n .GT. XMAX) XMAX = xm $ I F( X ( I) . LT. X MI N ) XMIN=X(I) 
IF( Y( n .GT.YMAXI YMAX=Ym i IF(Y ( I ) .LT. YMIN I YMIN=Y(I» 
1 CONTINUE * X(N+1)=XMIN « Y(N+1)=YMIN 

X(N+2)=.062E*(XMAX-XMIN) $ Y ( N + 2 » =. 0625* ( YM A X- YMI N ) 

CALL LINPLT(X,Y,N,K, -1,13, 3,0) 

CALL LINSEQI X ,Y,N,K,XLIM,YLIM,0,0,S I,0,0,0, 0,0) 

CALL LINS50(X,Y,N,K,XLIM,YLIM,1,0,SI,1,1,3,0,0) 

CALL LINSEQ(X,Y,N,K,XLIM,YLIM,8,0,SI,1,1,6,1,2) 

CALL CALPLT((l.,0.,<?oc) 


2 

FORMAT! 14) 

3 

FORMAT! 2F5 
END 

22 

1.0 

5.0 

2.0 

7.0 

2.2 

9.0 

2.5 

8.7 

3.0 

10. 0 

4.0 

11.0 

4.7 

12.C 

4.8 

11.7 

5.3 

12. 5 

6.5 

11.5 

7.2 

11.5 

7.5 

9.0 

S.O 

8.5 

8.5 

6.5 

9.0 

6.0 

9.5 

6.0 

9.7 

5.C 

10.0 

4.0 

11.0 

4.0 

11.2 

5.0 

11.5 

4.5 

12.0 

6.0 


Example 6 — Closed Curves 

Example 6 shows how CONSEQ can be used to draw closed curves. 

Result .- A four-point curve resembling a circle defined at (-1,0), (0,1), (1,0), and 
(0,-1) was drawn with Line-sequence parameters L = (0,0, 0,0,0). A four-point curve 
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resembling an ellipse defined at (-0.8,0), (0,0.5), (0.8,0), and (0,-0.5) was drawn with 
L = (1,1, 3, 0,0). These curves are presented below: 



Listing of program .- The program listing for Example 6 is as follows: 


DIMENSION X(23», Y(23),XLIM(2»fYLIM(2) 
♦ EXAMPLE 6, CLOSED CURVES 


ADJMIN=0. 

$ SCAFAC=. 

125 

CALL CALCOMP $ CALL 

LEROY $ CALL CALPLTCO 

SI=.16 $ 

K=1 



XLIMm=YLIM(l)=-1.2 

$ 

XLIM(^>=YLIM(2)=1.2 

N=4 




X(N+1)=Y(N+1 )=ADJMIN 

t 

X(N+2)=Y(N+2)=SCAFAC 

♦ FOUR POINT CIRCLE 



xm=o $ 

Ym=i, 



X(2)=l. $ 

Y(2)=0 



X(3»=0 i 

Y(3)=-l. 



X(4)=-l. $ 

Y(4)=0 




CALL CONSEQI X,y ,N,K,XLIM,YLIM,0,0,SIf 0f0f0,0,0) 
♦ FOUR POINT ELLIPSE 

X{ll=-,8 $ Y(1)=0 

X(2)=0 i Y(2)=.5 

X(3>=.8 $ Y(3)=0 

X(4)=0 i Y(4)=-.5 

CALL C0NSE0(X,Y,N,K,XLIM,YLIM,0,0,SI ,1,1 f3t OtO) 

CALL CALPLTCO. ,0,,999I 

END 
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Example 7 — Offset Feature 

Example 7 illustrates the use of the offset parameter to draw a heavy line curve 
such as would be used for an oral presentation. LINSEQ was called seven times with 
values of tlie offset parameter H differing by slightly less than the pen width. Note that 
for all calls except the first, the third parameter of LINSEQ was set to zero. 

Result .- The resulting figure plotted below is for the same input data points as 
Example 3: 



Listing of program .- The program listing for Example 7 is as follows: 

DIMENSION X(23),Y(23) ,XLIM(2),YLIM(2) 

* EXAMPLE 7. OFFSET FEATURE 

ADJMIN=0. $ SCAFAC=.l25 

CALL CALCOMP $ CALL LEROY $ CALL CALPLT (0. ,0. ,-3) 
SI=.16 $ K=1 $ N=2l 

X(N+1)=Y(N+1)=A0JMIN i X(N+2)=Y(N+2)=SCAFAC 
XLIM(1)=YLIM( 1)=-1.2 $ XLIM(2)=YLIM(2) = 1.2 

PI=3. 14159265358979 
no 1 1=1, N $ xn»=-l. + .l*(I-I) 

1 Y(I )=SIN(2.*PI*X(I) ) 

H=-.06 $ 00 2 11=1,7 

CALL LINSE0(X,Y,N,K,XLIM,YLIM,0,H,-SI,1, 1,6,1,2) 

H=H+.02 $ N=0 

2 CONTINUE 

CALL CALPLT(0.,0.,999) 

END 
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APPENDIX C 


SUBPROGRAM LISTINGS 

This appendix contains FORTRAN IV listings of subprograms LINSEQ and CONSEQ 
and of the subprograms called by LINSEQ and CONSEQ (except system routines CALPLT 
and SQRT). Annotation is minimal, so the descriptions in appendix A should be consulted 
when the listings are read. 
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Listing of Subroutine LINSEQ 


SUBROUTINE LINSEOl XX , YY , NN ,K , XL IM, YL I M, NSP , H, S I t L 1 , L2 , L3 ,L^ ,L5 ) 

* SUBROUTINE EQR PLOTTING SEQUENCE OF CURVED LINES USING LONG AND SHORT 
« DASHES 

DIMENSION XX( l),YY( 1) ,XLIM(1 »,YLIM( 1) 

COMMON/ A RCPAPVSI 101) / SC ALEX/ X( 101)/ SC ALEY/Y ( 10 1 ) 

COMMON/SLOPEX/OX( 101 ) /SLOPEY/OY ( 101 ) 

COMMON/C SP I SET/ JMAX , I FRMS, RMS 
JMAX=6 

IF( NN.EQ.O) GO TO 8 
N=NN 

kn-=k*n+i 

KNK=KN+K 

RX=1./XX(KNK) 

RY=1./YY(KNK) 

* COMPUTE SCALED FRAME LIMITS FOR LIMPLT 

X=( XLIM-XX( KN) )*RX 
X(2)=(XLIM( 2)-XX(KN) )*RX 
Y = ( YLIM-YY( KN) )*PY 
Y(2)=(YLIM(2)-YY(KN) ) *RY 
CALL LIMSET ( X, Y,0) 

* COMPUTE SCALED X,Y ARRAYS 

00 2 1=1, N 
K I=K*I 

X( I )=(XX(KI )-XX(KN) )*RX 
2 Y( I ) = (YY(KI )-YY(KN) )«RY 

* CALL PPETRP IF PRE-INTERPOLATION IS REQUESTED (SI.LT.O) 

IF( SI.LT.O)CALL PRETRP(N,SI) 

SS=SI 

* CALL ARCALC TO CCMPUTH APPROXIMATE ARC LENGTH Bt TWEEN DATA POINTS 

call ARCALC(N,0) 

IF(NSP.EQ.O) GO TO 6 

* SMOOTHING OF X = X(S) AMD Y = Y(S) IS PERFORMED BY CALLS TO SMOOTHC 

DO A 1=1, NSP 
0X= I 

CALL SMOOTHC (N,S,X,DX) 

DX= I 

A CALL SMOOTHC(N,S,Y,DX) 

CALL ARCALC(N,0) 

* CALLS TO SPISET INITIALIZE SPLINE DERIVATIVE ARRAYS DX AND DY FOR USE 
A BY SPIFIJN AND SPIDER 

6 CALL SPISET(N,S,X,DX) 

CALL SPISET (N,S,Y,OY) 

A CALL TO ARCSET INITIALIZES PEN POSITION 
8 CALL ARCSET(H,N,0) 

A CONFIGURE DESIRED SEQUENCE OF LONG AND SHORT DASHES WITH CALLS TO 
A ARCPLT WHICH EVALUATES INTERPOLATED CURVE AND CAUSES PEN MOTION 
IF( LI .EC.OGO TO la 
IFI LA.EQ.O) GO TO 20 
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* CALCULAT'^ REVISED POINT SPACING AND NUMBER OF CYCLES FOR INTEGER 

* NUMBER GF CYCLES PLUS TAIL OF LONG DASHES 

NIC = L2»'(L3 + L1 )+L4*(L5-*-Lll 

NIT = L2*( L3+LH-LI 

NIN=S(NJ /SS 

NCYC=N1N/NIC 

NIN=NCYC*NIC+MIT 

DS=S(N»/NIN 

DL1=L1*DS 

* general LINE-SFOUENCE CYCLE LOOP IS 14-LOOP 

no 14 ICYC=1,NCYC 
DO 12 12=1, L2 
CALL ARCPLT(OS,L3,2» 

12 CALL ARCPLTIOLl, 1,3 ) 

DO 1^ 14=1, L4 

CALL ARCPLT(0S,LS,2» 

14 CALL ARCPLT{DL1,1,3) 

* LOOP FOR TERMINATING TAIL OF LONG DASHES IS 16-LOOP 

DO 16 12=1, L2 
CALL ARCPLT ( DS,L3,2 ) 

16 CALL ARCPLT(0L1,1,3) 

RETURN 

* PLOT SOLID LINE 

13 NIN=S(N)/SS 
DS=S(N)/NIM 

CALL ARCPLT(DS,NIN,2) 

RETURN 

* PLOT LINE OF LONG DASHES ONLY 
20 NIC=L3+L1 

NIN=S(N)/SS 

NCYC=NIN/NIC 

NIN=NCYC*NIC-L1 

DS=S(N)/NIN 

DL1=L1*DS 

DO 22 ICYC=1,NCYC 

CALL ARCPLT(0S,L3,2» 

22 CALL ARCPLT(DL1,1,3) 

RETURN 

END SUBROUTINE LINSEO 
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Listing of Subroutine ARCPLT 


SURPOUTINE ARCPLT(DS,L, IPEN> 

♦ CURVE IS CALCULATED FROM SPLINE AND ARC OF LONG OP SHORT DASH IS 

♦ PLOTTfO PEN UP OR DOWN 

COMMON/ ARCPAR/S(ini ) /SCALE X/XI ID 1) / SC ALEV/ Y ( 1 0 1 ) 
COMMON/SLOPFX/nX( lOU /SLOPEY/DY( IDl ) 

IFIH.NE.OI GO TO A 
DO 2 1=1, L 
SP=SP+DS 

XP=5PIFUN(SP,N,S,X,OX) 

YP=SPIFIJN( SP,N,S,Y,OY) 

2 CALL LIMPLT(XP,YP,IPEN) 

RETURN 

A DO 6 1=1, L 
SP= SP+OS 

XP=SPIFUN( SP, N,S,X, OX ) 

XOOT= SPIDER ( SP,N,S, X,0X) 

YP=SPIFUN(SP,N,S,Y,DYI 
YDOT=SPIDER ( SP,N,S,Y,DY) 

HS=H/SOP.T( XOOT<'XOOT+VOOT*YDOT) 

XP= XP + HS*YD0T 
YP= YP-HS’f'XDOT 
6 CALL LIMPLT ( XP,YP, IPFNI 

return 

ENTRY ARCSET 

♦ ENTRY POINT ARCSET INITIALIZES H AND N, AND MOVES RAISED PEN TO 

♦ BEGINNING OF CURVE 

H=0 S 
N=L 
SP = 0 

call LIMPLT(X,Y,3I 

IF(H.EO.O) RETURN 

HS=H/SORT(OX*OX+OY*DYl 

XP=X+HS*DY 

YP= Y-HS*DX 

CALL LIMPLT(XP,YP,3) 

RETURN 

END SUBROUTINE ARCPLT 
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Listing of Subroutine PRETRP 


SUBROUTINE PRETRP(N,SI) 

♦ N X,Y POINTS ARE INTERPOLATED AS Y = Y(XJ TO GIVE Z^N-l POINTS IF X 

♦ IS MONQTONICALLY INCREASING. IF NOT, POINTS ARE INTERPOLATED AS 

♦ X=X(Y) IF Y IS MONOTONICALLY INCREASING, OTHERWISE RETURN. 

COMMON/ARCPAR/SdOl l/SCALEX/X( 101 )/SCALEY/Y( 101 1 
COMMON/SLOPEX/DXdOU /SLOPEY/OYdOl) 

SI=-SI 
00 1 1^2, N 

IF( X( I) .LE.X( I-in GO TO 4 

1 CONTINUE 
DO 2 1=1, N 
Sd ) = X( I ) 

2 DY{ n=Ym 

CALL SPISET(N,S,OY,OX,R.MSXI 
DO 3 1=2, N 
11=1+1-2 
12=11+1 
X(I2»=Sm 
Y(I2)=DY( I) 

x(m=.5*(s( i-ii+s( n I 

3 Y(Il) = SPIFUN(XdU,N,S,OY,OX» 

N=2’+N-l 

RETURN 

4 DO 5 1 = 2, N 

IF( Y( I) .LE.Y( I-D) RETURN 

5 CONTINUE 
DO 6 1 = 1, N 
Sd »=Y( n 

6 Dx( n=x( n 

CALL SPISET ( N,S,DX, OY,RNSY» 

DO 7 1=2, N 
11= I+I-2 
12=11+1 
Y( I2) = S( I ) 

Xd2)=DX( I ) 

Yd 1 ) = .5*(S( I-l»+S( I ) ) 

7 X( Il) = SPIFUN(Y(in ,N,S,OX,DYI 
N=2»N-l 

RETURN 

END SUBROUTINE PRETRP 


51 



*«• 


APPENDIX C — Continued 


Listing of Subroutine SPISET 

SUBROUTINE SP I SET ( N , X ,Y , 0 1 


♦ THIS SUBROUTINE COMPUTES DERIVATIVES OF A NATURAL CUBIC SPLINE FOR * 

♦ SUBSEQUENT USE BY FUNCTION SUBPROGRAMS SPIFUN AND SPIDER. * 

♦ STATEMENT 3 IMPLEMENTS EON IAA», PAGE 163 OF RALSTON + WILF, VOL 2 * 

♦ WITH NOTATION CHANGES NOTED BELOW, * 

♦ D(I-l) FOR U( II * 

♦ W FOR OMEQA = 4. * ( 2 SQR T ( 3 . >) ♦ 

♦ Y(I) FOR G( II * 

X( I ) FOR 8(11 * 

THE DO LOOPS IN THE SUBPROGRAM PERFORM THE FOLLOWING TASKS, * 

♦ 1 LOOP REPLACES X BY DIFFERENCES AND Y BY DIVIDED DIFFERENCES. ’S' 

♦ 2 LOOP REPLACES X BY SUBOIAGONALS (EQN (43) OF REF), Y BY 3* SECOND 

♦ DIVIDED DIFFERENCES, AND PUTS STARTING SECOND DERIVATIVES IN 0. * 

♦ 3 LOOP SOLVES ECN (41) OF REF BY OVERRELAXATION USING JMAX ♦ 

♦ ITERATIONS (DEFAULT VALUE OF JMAX IS 20). * 

♦ 4 LOOP UNDOES 2 LOOP. * 

♦ 5 LOOP REPLACES SECOND DERIVATIVES BY FIRST DERIVATIVES (AVERAGE OF ♦ 

♦ left and RIGHT). * 

♦ 7 LOOP COMPUTES RMS JUMP IN SECOND DERIVATIVES IF IFRMS=1. * 

♦ 9 LOOP UNDOES I LOOP. * 


DIMENSION X( 1 ) ,Y(1) ,D(1) 

COMMCN/CSPISET/JMAX, IFRMS.RMS 

DATA JMAX, IFPMS,R6,W/20,0, .166666656666667, 1.07179676972449/ 
0 ( 1 )= 0 ( N )=0 
00 1 11=2, N 
I=N+2-I I 

X( I )=X( I )-X( I-l) 

1 Y(I )=(Y( n-Y(I-l))/X( I) 

IF( N.EO.2) go TO 5 

DO 2 11=3, N 
I=N+3-I I 

D( I-l ) = 2.*( Y( I)-Y( I-l))/(Xm+X(I-D) 

Y( I )=1.5*D( I-l ) 

2 X( I )=.5*X( I-l )/(X( I )+X( I-l) I 
on 3 J=1,JMAX 

DO 3 1=3, N 

3 0( I-l )=W<‘( Y( I )-X( I ) ( I-2)-( .5-X( I) )*D( 1 ) )-(W-l. ) *D( I-l) 

DO 4 1=3, N 

x( I )=x( I-l ) ♦( .6/x( n-1. ) 

4 y ( I )=,333333333333333*Y( I )*(X(I)+X(I-1))+Y(I-1) 

5 SAVE=Y(2)-R6*X(2)*(2.*D(1)+D(2) ) 
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DO 6 1 = 2, N 
C=R6*X( I ) 

S1=Y( n-C<=( 2.*0( I-n+D( D) 

S2=Y( n+C*( D( I-l )+2 .*D( I » I 
D( I-l )=.5*( SAV=+S1) 

6 SAVE=S2 
D(N )=SAV6 

IF( IFRMS.NP.l) GO TO 8 
RMS=0 

IF(N.FQ.2) go to 8 
DO 7 1=3, N 
C=2.*D( I-l) 

7 RMS = RMS+{ (C+D( I-2)-3.*Y( I-l) )/X ( I-l ) + ( 0 ( I) +C-3 . *Y ( I ) ) /X( I) )**2 
RMS=2.*SQRT(RMS/N) 

8 DO 9 1=2, N 

Y(I ) = Y( I )*X( I )+Y( I-l ) 

9 X(I )=X( I)+X( I-l) 

RETURN 
FNO 


SUBROUTINE SPISET 


appendix C — Continued 

Listing of Function Subprogram SPIFUN 

FUNCTIOM SPIFUNi XP,N,X,Y,0> 

DIMENSION X( 1 ) ,Y ( 1) ,C( 1 ) 

EVALUATES A NATURAL CUBIC SPLINE USING SLOPE ARRAY 0 CALCULATED BY 
SPISET AND USING THE INPUT DATA ARRAYS X AND Y. 

IF( XP.LE.XIl ) » GC TO 6 
DO 2 1^2, N 

IF(XP.LT.X( I n GO TO 4 
2 CONTINUE 

SPIFUN=Y(N)+D(N)*( XP-X(N)» 

RETURN 

4 ci=i./(x( n-x(i-in 

C2=X( I )-XP 
C3=XP-X( I-l ) 

C4=C2*C1 

C5=C3*C1 

SPIFUN=C5*C5<'( ( l.+2.*C4 )*Y(I »-C2*D( I ) ) 

, +C4«C4*((l.+2.<‘C5)*Y( I-U+C3*0( I-l) ) 

RETURN 

6 SPIFUN=Y (1 |-D(1 ) *( X( 1 )-XP) 

RETURN 

END REAL FUNCTION SPIFUN 


Listing of Function Subprogram SPIDER 


FUNCTION SPIDER (XP,N,X,Y,0) 

DIMENSION X (1 ) ,Y(L) ,D( 1 ) 

EVALUATES THE FIRST DERIVATIVE OF A NATURAL CUBIC SPLINE USING SLOPE 
ARRAY D CALCULATED BY SPISET AND USING THE INPUT DATA ARRAYS X AND Y 
IF( XP.LE.XI I n GO TO 6 
DO 2 1=2, N 

IF( XP.LT .XU n GO TO 4 
2 CONTINUE 
SPIDER=D(N) 

RETURN 

4 ci=i./( x( n-x( i-in 
C2=X( I )-XP 
C3=XP-X( I-l) 

C4=2.*C2-C3 

C5=2.*C3-C2 

S P I DE R= C 1 *C 1 « ( C 3* ( 2 . * ( 1 . «-C 1 *C4 ) * 

, Y( I )-C4#0( I) )-C2*( 2.*( l.^CI*C5)*Y{ I-1H-C5*D( I-IH) 

RETURN 

6 SPIOER=Dm 
RETURN 
END 


REAL FUNCTION SPIDER 



APPENDIX C — Continued 


Listing of Subroutine SMOOTHC 


SUBROUTINE SMOOTHC ( NtX ,Y,T ) 

* THE Y ARRAY IS SMOOTHED BY A LOCAL FIVE POINT LEAST SQUARES CUBIC 

♦ WEIGHTED BY W 

DIMENSION Xm,Y(l»,T(ll - 

IF(N.LT,5)RETURN 
S=(T*(X(N»-Xm )/N)*#2 
DO A L=ltN 

K=MIN0(N-4tMAX0( l.L-21 ) 

K4=K+4 

00 1 1 = 1,20 

1 T(I) = 0. 

DO 3 M=K,K4 

W=l,/( S+(X(L )-X(Mn*»2 ) 

R=l. 

DO 3 1=1,4 
RR=1. 

DO 2 J=l,4 

T( I+J*4-4)=T ( l+J*4-4)+R*RR*W 

2 RR=RR*X(M> 

T( I+16)=T(I*16)+R*Y(MJ*M 

3 R=R*X(M) 

CALL CHLSKYS(T,4,T(17» ,1,41 
M=1+M00(L-1,5» 

IF(L.GT,5) Y(L-5)=T(M+20» 

T(M+20I=0, 

R = l. 

DO 4 J=l,4 

T(M+20)=T(M+20»+R*T( J+16) 

4 R=R*X(U 
DO 5 L=1,5 

ML=l + MOO(M+L-l , 5 ) 

5 Y(N-5+L)=TIML+20) 

RETURN 

END SUBROUTINE SMOOTHC 
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Listing of Subroutine LIMPLT 

SUBROUTINE L IMPLT ( X , Y , I PEN » 

♦ CAUSES PEN MOTION IF (X,Y» IS IN PREDECLARED FRAME 

DIMENSION Xm,Y(l) 

IPLM=IPEN 

IF(X.GT.XX) GO TO 1 
IF(X.LT.XN) GO TO 1 
IF(Y.GT.YX) GO TO 1 
IFIY.LT. YNI GO TO 1 
IF< IFLAG.EQ. 1) IPLM=3 
CALL CALPLT(X,Y, IPLM) 

IFLAG=0 

RETURN 

1 IFLAG=1 
RETURN 

ENTRY LIMSET 

♦ ENTRY POINT LIMSET IS USED TO PASS FRAME LIMITS (INCHES) 

IFLAG=IPEN 

XN=X 

YN=Y 

XX=X(2) 

YX=Y( 21 
RETURN 

END SUBROUTINE LIMPLT 

Listing of Subroutine CHLSKYS 

SUBROUTINE CHLSKYS ( A ,N f B,M ,NX) 

DIMENSION A(NX,1) ,e(NXfl) 

♦ CHOLESKY DECOMPOSITION IS USED TO SOLVE THE MATRIX EQUATION AX = B, 

♦ WHERE THE COEFFICIENT MATRIX, A, IS SYMMETRIC. ON OUTPUT X IS 

♦ STORED IN B, 

IF(N.EQ.l) GO TO 6 
DO 2 1=2, N 
11 = 1-1 
DO 2 J=I ,N 
DO 2 L=1,I1 

2 All , J)=A( I , J )-A( L, n*A (L,J)/A(L,L) 

DO 5 K = 1 ,M 

DO 3 1=2, N 
11 = 1-1 
DO 3 L = 1 ,11 

3 B(I,K)=R(I,K)-A(L,I)*8(L,K)/A(L,L) 

DO 4 1=2, N 

11 = 1-1 
DO 4 L=1 , 1 1 
NI=N-I1 
NL=N+1-L 

4 B(NI,K)=B(NI ,K)-A(NI,NL)*B(NL,K)/A(NL,NL) 

DO 5 1=1 ,N 

5 B( I ,K)=B( I ,K )/A( I, I ) 

RETURN 

6 A(l,l)=l./A(l,l) 

DO 7 L = 1 ,M 

7 B(1,L)=A(1,1)*B(1,L) 

RETURN 

END SUBROUTINE CHLSKYS 
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Listing of Subroutine ARCALC 


SUBROUTINE ARCALCC N t IFLAGJ 

I* ARCALC CALCULATES THE S(l) ARRAY OF APPROXIMATE ARC LENGTHS 
COMMCN/ARCPAR/S( 101 )/SCALEX/X{ 101 )/ SCALE Y/Y( 101 » 

R24=.0A1 66666666666667 

sa )=o 

S(2 ) = B1=SQRT( ( X( 2)-X( 1 ) ) *<= 2 + ( Y ( 2 J-Y ( 1 >>*<‘2 ) 

IF ( N. EQ. 2> RETURN 

N1=N-1 

HS=H1=0 

IF( IFLAG.EQ.O) GO TO 1 

B2 = SQRT( (xm-X(Nin**2+(Y(l)-Y(Nl» 

C1 = SQRT( (X(2)-X(N1) )**2+(Y(2)-Y{Nin**2) 

Al = (X(U-X(Nl J »*( Y( 2I-Y(N1 n-(X( 21-XINl n«( Y( 1 l-YINl I I 
HS=Hl=Al/( B2*B1*C1 » 

1 DO 2 1=2, N1 
IP=I+1 
IN=I-1 

B=SQRT((X(IP)-X( n )**2 + (Y( IP)-Y(m<'*2» 

C=SORT((XHP J-X( INn*«‘2 + (Y( IP)-y(INn**2) 

A=( X( I )-X( IN n*( Y( IP|-Y( IN) )-(X( IP>-X( IN)) *( Y( n-Y( IN) I 
H=A/(B1*B*C) 

S( I )=S( I-D + Bl* ( 1. +R2A*( (Hl+HM-Bl )#*2 ) 

Bl = B 
H1 = H 

2 CONTINUE 

S(N)=S(N-1)+B1*( l.+R24*( (H1+HS)*B1)*«2) 

RETURN 

ENO SUBROUTINE ARCALC 
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Listing of Subroutine GONSEQ 


subroutine CONSeQ( XXtYY, NN,K,XLIM,YLIM,NSP, H, SI ,Ll,L2,L3,L^tL5 » 

* SUBROUTINE FOR PLOTTING SEQUENCE OF SMOOTH CLOSED CURVES (CONTOURS) 

* USING LONG AND SHORT DASHES 

DIMENSION XX( 1), YY( U ,XL IM(H,YLIM( 1) 

COMMON/ARCP AR/S( 101 )/ SC A LE X/X(l 01 )/ SC ALEY/Y { 1 0 1 ) 

COM MON/ SLOPE X /OX ( 101 ) /SLOPEY/OY( 101 ) 

COMMON/CSPI SET/JMAX, I FRMS.RMS 
JMAX=6 

IF( NN.EQ.O)GO TO 3 

N = NN-H 

SS=SI 

KN=K*NN+1 

KNK = KN+K 

RX=1./XX(KNK) 

RY=1 . /YY( KNK ) 

* COMPUTE scaled FRAME LIMITS FOR LIMPLT 

X=( XLIM-XX< KN) )«RX 
X(2)=(XLIM(2)-XX(KN) ) *RX 
Y=( YLIM-YY(KN) )*RY 
Y (2) = (YLIm{ 2)-YY(KNn*RY 
CALL HMSET(X,Y,0) 

* COMPUTE SCALED X,Y ARRAYS 

DO 2 1=1, NN 
KI=K*I 

X( I ) = (XX(KI )-XX(KN) I'i'RX 
2 Y(I )=(YY(Kl)-YY(KN) )*RY 
X(N)=X(1) 

Y(N )=Y( 1 ) 

* CALL ARCALC TO COMPUTE APPROXIMATE ARC LENGTH BETWEEN DATA POINTS 

CALL ARCALC(N,1) 

TF(NSP.EO.O) GO TO 6 

* SMOOTHING OF X = X(S) AND Y = YIS) IS PERFORMED BY CALLS TO SMOOTHP 

DO A I=1,NSP 
0X= I 

CALL SMCOTHP(N,S,X,DX) 

DX= I 

A CALL SMOOTHP(N,S,Y,DX,DY) 

CALL ARCALC (N,l) 

* CALLS TO SPISETP INITIALIZE SPLINE DERIVATIVE ARRAYS DX AND OY FOR US 

* BY SPIFUNP AND SPIDERP 
6 CALL SPISETP(N,S,X,DX) 

CALL SPISETP(N,S,Y,DY) 

* CALL TO ARCSETP INITIALIZES PEN POSITION 

8 CALL ARCSETP(H,N,0) 

* CONFIGURE DESIRED SEQUENCE OF LONG AND SHORT DASHES WITH CALLS TO 

ARCPLTP WHICH EVALUATES INTERPOLATED CURVE AND CAUSES PEN MOTION 
IF( Ll.EO.OIGO TO 13 
IF(LA.E0,0)G0 TO 20 
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* CALCULATE REVISED POINT SPACING AND NUMBER OF CYCLES FOR INTEGER 

* NUMBER OF CYCLES 

NIC=L2*(L3+L1I+LA*(L5+L1J 

NIN=S(N) /SS 

NCYC=NIN/NIC 

NIN=NCYC*NIC+L1 

OS=S(Nf /NIN 

DLl=Ll*OS 

» GENERAL LINE-SEQUENCE CYCLE LOOP IS 14-LOOP 
DO ICYC=1,NCYC 
DO 12 12=1, L2 
CALL ARCPLTP(DS,L3,2I 
12 CALL ARCPLTP(DL1,1,3) 

DO 14 14=1, L4 

CALL ARCPLTP(OS,L5,2) 

14 CALL ARCPLTP(DL1,1,3) 

RETURN 

* PLOT SOLID LINE 
18 NIN=S(N)/SS 

OS=S(NI /NIN 

CALL ARCPLTP(DS,NIN,2 ) 

RETURN 

* PLOT LINE OF LONG DASHES ONLY 
20 NIC=L3+L1 

NIN=S(N» /SS 
NCYC=NIN/NIC 
NIN=NCYC*NIC 
DS=S(N)/NIN 
0L1=L1*DS 
00 22 ICYC=l,NCYC 
CALL ARCPLTP(0S,L3, 21 
22 CALL ARCPLTP(0L1,1,3) 

RETURN 

END SUBROUTINE CONScO 
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Listing of Subroutine ARCPLTP 


SUBROUTINE ARCPLTP ( DS t L , IPEN ) 

♦ CURVE IS CALCULATED FROM SPLINE AND ARC OF LONG OR SHORT DASH IS 

♦ PLOTTED PEN UP OR DOWN 

COMMON/ ARCPAR/S ( lOU /SC ALEX/XI 1 01 »/ SC ALE Y/ Y ( 1 0 1 ) 
COMMCN/SLOPEX/DXdOl l/SL0PEY/0Y(101 ) 

IF(H.NE.O) GO TO 4 
DO 2 1=1 ,L 
SP=SP+DS 

XP = SPIFUNP(SPtN,S,X,DX ) 

YP=SPIFUNP(SP,N,StY,CY» 

2 CALL LIMPLTI XP,YP, IPEN ) 

RETURN 

4 DO 6 1=1, L 
SP^SP+DS 

XP=SPIFUNP (SP,N, S, X,OX I 
XDOT=SPIDERP(SP,N,S ,X, CX 1 
YP=SPIFUNP(SP,N,S,Y,CY) 

YDOT=SPIOERP( SP,N, S,Y,DY) 

HS = H/SQRT(XOOT*XDOT + YDC.T*YDOT) 

XP=XP+HS*YDOT 
YP=YP-HS*XDOT 
6 CALL LIMPLTI XP, YP, IPEN ) 

RETURN 

ENTRY ARCSETP 

* ENTRY POINT ARCSETP INITIALIZES H AND N, AND MOVES RAISED PEN TO 

* BEGINNING OF CURVE 

H=DS 
N=L 
SP = 0 

CALL LIMPLTI X,Y ,3 ) 

IFIH.EQ.O) RETURN 
HS=H/SQRT I OX*DX+CY*OY ) 

XP=X+HS*DY 

YP=Y-HS*DX 

CALL LIMPLTIXP,YP,3) 

RETURN 

END SUBROUTINE ARCPLTP 
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Listing of Subroutine SPISETP 


SUBROUTINE SPISeTP(N,X,Y,D) 

** >![***:«t** ijt :{£:<( 4c 4= A 4t 4t* ’ll ♦**#* "tt *♦♦*♦##<'** * 


* THIS SUBPQUTINF COMPUTES DERIVATIVES OF A PERIODIC CUBIC SPLINE FOR 

* SUBSEQUENT USE BY FUNCTION SUBPROGRAMS SPIFIJNP AND SPIOERP, 

* THE SUBROUTINE SETS YIN) TO Y(l). THE PERIOD IS X(N)-XIl). 

* statement 4 IMPLEMENTS EON 1^4), PAGE 143 OF RALSTON + WILF, VOL 2 

* WITH NOTATION CHANGES NOTED BELOW, 

* DII-1) FOR UI I) 

* W FOR OMEQA = 4.4(2.-S0RT{3. )) 

* Yd) FOR G( I ) 

* B FOR 8(1) 

* THE DO LOOPS IN THE SUBPROGRAM PERFORM THE FOLLOWING TASKS, 

* 2 LOOP PUTS STARTING SFCQNO DERIVATIVES IN D. 

* 5 LOOP SOLVES EON (41) OF REF BY OVERREL AXAT ION USING UMAX 

* ITERATIONS (DEFAULT VALUE OF UMAX IS 20. 

* 3 LOOP REPLACES SECOND DERIVATIVES BY FIRST DERIVATIVES (AVERAGE OF 

* LEFT AND RIGHT). 

* 10 LOOP COMPUTES RMS JUMP IN SECOND DERIVATIVES IR 1FRMS=1. 


A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 


AA A A AAA AA A AAA A A AA A A A A A A A A A AAA A A A A A AAA AAA A A A A A AAA A AAA A A A AA AA A A A A A AAA A AA A A A A A 


dimension X(1),Y(1) ,D(1) _ 

COMMON/CSPI SET /UMAX, I FR MS, RMS 

DATA UMAX, I FRMS , R6 , W/20 , 0 ,. 164666666666667 , 1. 0T179676972449/ 
Y(N)=Y( 1 ) 

P=X(N)-X( 1) 

IF(N.EQ,2) GO TO 12 
N1=N-1 

0(1 ) = 2.*( (Y(2 )-Y( in/(X(2)-X( 1) )-(Y(U-Y(N-l) ) / ( X(1 )-X(N-l) + P) )/ 

, (X(2)-X(N-1 )+P) 

DO 2 1=2, N1 

2 D( I 1 = 2. *( ( Y{ lAl )-Y( I) ) / (X( I+U-X( I) )-( Y( I)-Y( I-l) )/(X( I)-X( I-l) ) )/ 
, (Xd+1)-X( I-l) ) 

DIN ) = D( 1 ) 

DO 6 J=1,JMAX 

B=.5A(X( I )-X(N-l)+P)/(X(2)-X(N-l)+P) 

D(l )=W*(3.*((Y(2)-Y(l n/(X(2l-X(lM-(Y(l)-Y(N-l) )/(X(l)-X(N-ll+Pn 
, /( X(2)-X( N-1) AP)-B*D (N-1 )-( .5-B)#0(2) )-(W-l. )*D( 1) 

DO 4 1=2, N1 

B=.5*(x( n-x(i-i))/(x(i+i)-x(i-D) 

4 D( I )=W*(3.*( (Y( I + l)-Y( I ))/(X(lAl)-xm)-(Y(I)-Y(I-l))/(Xm-X(I-l) 
, ) ) / (X( l+l )-X(I-in-B*0( I- l)-( .5-B)*D( l + l) )-( W-1. )A0( I ) 

6 D(N)=D(1) 
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SAVe^(Y(2)-Y(l))/(X(2)-X( 1) X( 2>-X( U ) *( 2 . ( 1» +0 ( 2 ») 

DO a 1=1, N1 

YDIF=(Y(I+l)-Y(I)(/(X(I + l)-X(in 
X0IF=R6>!'(X( I + l)-X( I ) ) 

Sl=YDIF-XDIF«(2.<=0( I + l) ) 

S2=Y0!F + XDIF*{D( I > + 2. *D( I + IU 
0(1 ) = .F*( SAV£: + S1 ) 

3 SAVE=S2 
D(NJ=0(1) 

IF( IFRMS.NE.l ) RFTURM 

si=i./( x( n-x(N-i > + p) 

S2=l./(X(2I-X(in 

RMS=( Sl*( 2. *D( 1 )+U(N- U-3. *Sl*( Y( 1)-Y( N-l> ) )^S2*(D( 2)+2.* 

, Dl 1 )-3.*S2«( Y(2)-Y(l n ) »>^*2 
DO 10 1=2, N1 
si=i./(x( n-x(i-D) 

$2= 1 ./( X( I + l »-X( I ) » 

10 RMS = PMS+(S1*( 2.*0( n+0( I- 1 )-3.*Sl*( Y(I )-Y( I-l ) )■) +S2*(0( I^lH-2. * 
, D( I )-3.*S2*( Y(!+l )-Y( I n ) )**2 
P*^S = 2.*SQRT (PMS/(N-1. ) » 

RFTUPM 

120(1 )=D( 2 )=0 
RETURN 

END SUBROUTINE SPISETP 
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Listing of Function Subprogram SPIFUNP 

FUNCTION SPIFUNP(XP,N,X,Y,D) 

♦ EVALUATES A PERIODIC CUBIC SPLINE USING SLOPE ARRAY D CALCULATED BY 

♦ SPISETP AND USING THE INPUT DATA ARRAYS X AND Y 

DIMENSION xm,Y(l) ,0(11 
XX»XP 

P=X{NJ-X(1I 

2 IF(XX.GE.X(1 n GO TO 4 
XXsXX+P 
GO TO 2 

4 IF(XX.LT.X(N ) ) GO TO 6 
XXsXX-P 
GO TO 4 
6 DO 8 I=2,N 

IF( XX,LT.X( I n GO TO 10 
8 CONTINUE 

10 C1=1./(X( I »-X(I-l) ) 

C2=X( I )-XX 
C3=XX-X( I-l) 

C4=C2*C1 

C5=C3*C1 

SPIFUNP=C5*C5*( ( 1. +2.*C4)*Y( I )-C2*D( I ) ) 

, +C4*C4*((1.+2.*C5I*Y( I-1)+C3*D( I-l) ) 

RETURN 

END REAL FUNCTION SPIFUNP 


Listing of Function Subprogram SPIDERP 


FUNCTION SPIDERP(XP,N,X,Y,0) 

* EVALUATES THE FIRST DERIVATIVE OF A PERIODIC CUBIC SPLINE USING SLOPE 

♦ ARRAY D CALCULATED BY SPISETP AND USING THE INPUT DATA ARRAYS X AND Y 

DIMENSION X(l ),Y(1) ,D(1) 

XX=XP 

p=X (N)-X( 1 ) 

2 IF( XX.GE.XI 1) ) GO TO 4 
XX=XX+P 
GO TO 2 

4 IF( XX.LT.X( N) ) GO TO 6 
XX=XX-P 
GO TO 4 
6 DO 8 1 = 2, N 

IF( XX.LT.X( I) ) GO TO 10 
8 CONTINUE 

10 Cl=l./( X( I )-X( I-l) ) 

C2=X( I )-XX 
C3=XX-X( I-l ) 

C4=2.*C2-C3 

C5=2.t'C3-C2 

SPIDERP = C1*C1*(C3*( 2.*( l.+Cl*C4)*' 

, Y( I )-C4<'D( I ) )-C2*(2.*( l.+Cl«C5)*Y( I- 1 ) +C5 *D( I - 1 ) ) ) 

RETURN 

END REAL FUNCTION SPIDERP 
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Listing of Subroutine SMOOTHP 

SUBROUTINE SMOOTHP ( N t X ,Y»T I 
DIMENSION xm , Y(1 » ,T( U ,SAVE(3) 

* THE Y ARRAY IS SMOOTHED BY A LOCAL FIVE POINT LEAST SQUARES CUBIC 

♦ WEIGHTED BY W (Y IS A PERIODIC FUNCTION OF X) 

IFIN.LT.SJRETURN 
P=X(N)-X(1 ) 

SAVE(2)=Y(2I 

SAV£(3>=Y(3) 

S=(T*(X(N)-X(in/N)**2 . 

DO 6 L=1 f N 
DO 1 1=1,20 

1 T(I»=0, 

DO 5 K=l,5 

M=L+K-3 

XM=XIM) 

YM=Y(M) 

IFIM.GT.OJ GO TO 2 

MM=M+N-1 

YM=Y(MM) 

XM=X(MM)-P 

2 IF(M.LE.N) GO TO 3 
MM=M+1-N 
YM=SAVE( MM) 

XM=X(MM) +P 

3 W=1./(S+{XIL)-XM)**2) 

R=l, 

DO 5 1=1 ,4 
RR=1. 

DO 4 J=l,4 

T( I+J*4-4) = T( I+J*4-4)+P*RR*W 

4 RR=RR*XM 

T( I +16I=T{ I +16) +R’«=YM*W 

5 R=R*XM 

CALL CHLSKYS(T,4,T( 17) ,1,4) 

M=1+M0D(L-1, 5) 

IF<L.GT.5) YIL-5 >=T(M+20) 

T(M+20)=0, 

R=l. 

DO 6 J = 1 ,4 

T( M+20)=T(M+20)+R*T( J+16) 

6 R=R*X(L) 

DO 7 L=lf5 
ML=1+M0DIM+L-1,5) 

7 Y(N-5+L)=T(ML+20) 

RETURN 

END SUBROUTINE SMOOTHP 
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edge of phenomena in the atmosphere and space. The Administration 
shall provide for the widest practicable and appropriate dissemination 
of information concerning its activities and the results thereof." 

— National Aeronautics and Space Act of 1958 


NASA SCIENTIFIC AND TECHNICAL PUBLICATIONS 


TECHNICAL REPORTS: Scientific and 
technical information considered important, 
complete, and a lasting contribution to existing 
knowledge. 

TECHNICAL NOTES: Information less broad 
in scope but nevertheless of importance as a 
contribution to existing knowledge. 

TECHNICAL MEMORANDUMS: 
Information receiving limited distribution 
because of preliminary data, security classifica- 
tion, or other reasons. 

CONTRACTOR REPORTS: Scientific and 
technical information generated under a NASA 
contract or grant and considered an important 
contribution to existing knowledge. 


TECHNICAL TRANSLATIONS: Information 
published in a foreign language considered 
to merit NASA distribution in English. 

SPECIAL PUBLICATIONS: Information 
derived from or of value to NASA activities. 
Publications include conference proceedings, 
monographs, data compilations, handbooks, 
sourcebooks, and special bibliographies. 

TECHNOLOGY UTILIZATION 
PUBLICATIONS: Information on technology 
used by NASA that may be of particular 
interest in commercial and other non-aerospace 
applications. Publications include Tech Briefs, 
Technology Utilization Reports and 
Technology Surveys. 


Details on the availability of these publications may be obtained from: 

SCIENTIFIC AND TECHNICAL INFORMATION OFFICE 

NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 

Wathinglon, O.C. 20546 




